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ABSTRACT 


The  problem  of  imaging  through  atmospheric  turbulence  has  plagued  astronomers  for 
centuries.  Only  in  the  last  20  years  has  the  lost  resolution  been  restored  through  innova¬ 
tive  hardware  and  algorithms.  This  report  discusses  the  new  optical  imaging  technologies 
of  adaptive  optics  and  speckle  imaging  in  the  space  surveillance  context.  The  techniques 
are  compared  and  contrasted,  and  new  imaging  methods  are  proposed  which  combine  ele¬ 
ments  of  both.  New  methods  are  suggested  for  optical  data  evaluation. 
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L  INTRODUCTION 


Optical  study  of  space  objects  has  faced  atmospheric  limitations  to  resolution  for  centuries.  One 
finds  that  increasing  a  lens  diameter  beyond  about  10  cm  is  ineffective  in  increasing  telescope  resolution. 
Light-gathering  ability  increases,  but  atmospheric  turbulence  prevents  a  corresponding  increase  in 
resolution. 

In  the  last  15  years,  two  new  and  exciting  approaches  to  imaging  through  turbulence,  adaptive 
optics  and  speckle  imaging  have  proven  able  to  greatly  reduce  the  degrading  effects  of  turbulence  in 
some  cases.  This  offers  hope  that  ground-based  optical  sensors  may  become  contributors  to  the  high- 
resolution  imaging  of  space  objects,  an  area  now  dominated  by  radar.  There  may  also  be  opportunities 
for  data  fusion,  combining  information  from  optical  and  radar  sensors. 

One  current  sensor  of  interest  is  the  AMOS/MOTIF  (ARPA  Maui  Optical  Station/Maui  Optical 
Tracking  and  Identification  Facility)  Compensated  Imaging  System  (CIS).  This  adaptive  optics  device 
has  been  able  to  provide  high-resolution  satellite  images  in  some  cases. 

The  adaptive  optics  approach  requires  complicated  mechanical  equipment  in  a  feedback  control 
arrangement.  There  is  extensive  real-time  processing  required,  but  the  algorithms  are  relatively  straight¬ 
forward.  The  speckle  approach  requires  simpler  hardware,  but  utilizes  extensive  non-real-time  process¬ 
ing  with  innovative  and  complicated  algorithms. 

These  two  new  technologies  have  a  lot  in  common,  but  the  literature  is  fairly  disjoint.  This  report 
will  attempt  to  explain  the  relationships  between  them,  as  well  as  discuss  the  special  problems  of  satel¬ 
lite  surveillance.  There  is  some  synergy  between  these  areas  which  can  benefit  each  of  them. 

One  question  addressed  by  this  report  is  whether  or  not  the  algorithms  recently  developed  for 
speckle  processing  may  be  used  to  post-process  adaptive  optics  data,  increasing  the  effective  resolution. 
One  obviously  cannot  combine  these  techniques  in  the  other  order,  i.  e.,  once  data  have  been  collected 
for  speckle  processing  it  is  too  late  to  apply  adaptive  optics.  It  is  possible,  however,  that  the  adaptive 
optics  approach  may  suggest  new  or  improved  speckle  processing  algorithms. 

This  report  is  written  from  the  point  of  view  of  signal  processing  and  systems  engineering.  The 
required  background  for  this  report  includes  random  process  theory,  basic  digital  signal  processing,  and 
basic  physics. 

To  gain  a  complete  understanding  of  these  techniques  requires  background  in  physics,  optics, 
astronomy,  control  theory,  signal  processing,  mathematical  optimization,  and  computer  science,  among 
other  subjects.  Sections  2-6  of  this  report  will  present  a  fair  amount  of  background  material  on  optics, 
turbulence,  adaptive  optics,  and  speckle  processing  to  make  the  analysis  accessible  to  those  not  trained 
in  optics.  The  goal  is  to  provide  sufficient  detail  that  the  physics  of  the  situation  will  become  clear, 
quantitative  calculations  can  be  made,  and  the  reader  can  consult  the  literature  for  details  and  new  devel¬ 
opments.  Derivations  of  well-known  results  will  be  left  to  the  references. 

Section  7  will  attempt  to  present  a  more  unified  analysis  of  the  new  optical  technologies,  and  dis¬ 
cuss  satellite  imaging  problems  in  more  detail.  The  new  results  are  primarily  confined  to  this  section. 
These  include  an  analysis  of  adaptive  optics  errors  and  a  proposed  method  for  speckle  data  evaluation. 
The  reader  who  is  well  acquainted  with  speckle  processing  and  adaptive  optics  is  encouraged  to  start 
with  Section  7,  returning  to  previous  sections  only  as  necessary. 
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2.  OVERVIEW 


2.1  OPTICS 

In  radar  and  communications,  we  often  deal  with  coherent,  narrowband  systems.  We  can  envision 
waves  propagating  through  space,  and  directly  detect  and  process  them  coherently.  While  antennas  are 
often  two-dimensional,  the  signals  are  quickly  beamformed  into  individual  time  functions. 

In  optics,  we  often  deal  with  thermal  light,  e.  g.,  light  emitted  from  excited  atoms  in  some  random 
configuration.  The  resulting  signal  is  a  sample  from  a  random  process  which  is  wideband  and  non¬ 
coherent.  We  also  are  generally  interested  in  image  formation,  so  rather  than  time  functions  we  have 
functions  of  two  or  three  spatial  dimensions  plus  time. 

The  best  way  to  handle  such  signals  is  to  use  what  is  known  as  the  analytic  signal  representation. 
This  method  is  occasionally  used  in  the  study  of  modulation  and  demodulation  in  communications,  but 
is  more  prevalent  in  statistical  optics.  Our  goal  is  to  discuss  the  propagation  of  mutual  coherence  func¬ 
tions  and  the  important  Van  Cittert-Zemike  Theorem. 

2.2  PHYSICAL  MODEL 

The  physical  situations  of  interest  in  this  report  may  be  described  as  follows.  In  the  first  case,  an 
astronomical  object  emits  thermal  light.  In  the  second  case,  radiation  from  a  source  of  thermal  light 
(such  as  the  sun)  is  reflected  from  a  satellite.  In  both  cases,  the  light  travels  through  the  atmosphere  and 
is  modified  by  turbulent  layers  before  entering  the  aperture  of  a  telescope.  The  thermal  light  is  best 
modeled  as  a  random  process,  and  the  atmosphere  likewise  fluctuates  randomly. 

The  limiting  resolution  which  a  telescope  can  provide  is  approximately 


where  D  is  the  telescope  diameter  and  X  is  the  wavelength.  The  effects  of  turbulence  can  be  included  in 
this  expression  by  replacing  the  telescope  diameter  D  by  an  effective  D.  In  areas  of  good  astronomical 
“seeing”  a  long  film  exposure  gives  an  effective  D  of  about  10  cm.  For  example,  with  X  =  500  nm,  the 
resolvable  angle  0  is  510“^  rad,  or  about  1”  (second  of  arc).  The  light  collected  is  proportional  to  the 
area  of  the  telescope,  but  the  resolution  does  not  improve  as  D  increases  beyond  10  cm. 

Another  way  of  looking  at  this  is  in  the  frequency  domain.  If  we  take  the  two-dimensional  Fourier 
transform  of  such  an  image,  spatial  frequencies  greater  than  about  one  cycle/arc  sec  in  either  direction 
will  be  small. 

23  ADAPTIVE  OPTICS 

Adaptive  optics  has  been  discussed  and  even  implemented  as  long  ago  as  the  1950s.  The  idea  is 
not  new,  but  the  bandwidths  of  the  original  systems  were  in  the  region  of  several  Hz  instead  of  the  sev¬ 
eral  hundred  or  thousand  Hz  required  to  correct  atmospheric  turbulence.  Hardy’s  useful  survey'  dis¬ 
cusses  these  early  efforts.  The  first  reported  high-bandwidth  adaptive  optics  system  was  produced^  in 
1973.  The  monolithic  piezoelectric  mirror  (MPM)  is  the  device  which  mechanically  deforms  at  high 
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rates  in  response  to  the  controlling  voltage.  Interferometers  and  an  analog  computer  are  the  other  key 
components  in  this  system. 

Atmospheric  turbulence  will  be  described  in  greater  detail  below,  but  during  periods  of  good  see¬ 
ing,  the  major  effect  of  turbulence  is  phase  shifting  due  to  path  length  variations;  the  amplitudes  of  the 
light  are  not  greatly  affected.  In  theory,  mechanically  deforming  a  mirror  in  the  light  path  can  perfectly 
correct  the  phase.  Since  the  residual  amplitude  error  is  much  less,  we  approach  the  diffraction  limit  in 
resolution. 

2.4  SPECKLE  IMAGING 

In  1970  Labeyrie  published  a  landmark  paper^  ushering  in  the  era  of  speckle  processing.  Notation- 
ally,  let  /  (.r,y)  be  the  image,  and  F(u,v)  its  Fourier  transform.  Labeyrie’s  insight  was  that  although  one 
could  not  achieve  the  diffraction  limit  by  time  averaging  /  (x,y)  (long  film  exposure)  things  improved 
greatly  when  the  quantity  |  F  (u,v)  \  ^  was  time  averaged. 

If  we  average  the  spatial  frequencies  F  {u,v),  we  are  averaging  many  complex  quantities  with  zero 
mean  so  our  result  tends  toward  zero.  (Note  that  averaging  /  {x,y)  is  equivalent  to  averaging  F  {u,v),  by 
linearity.)  On  the  other  hand,  if  we  average  the  squared  magnitudes  of  the  spatial  frequencies  (the  power 
spectrum)  we  are  averaging  positive  quantities  so  the  resultant  does  not  tend  toward  zero. 

From  F(u,v),  one  could  retrieve  the  desired  high-resolution  image  using  an  inverse  Fourier  trans¬ 
form.  However,  only  I  F(m,v)  |  the  power  spectrum,  is  available.  The  inverse  transform  of 
I  F(m,v)  I  ^  is  the  autocorrelation  function  of  the  image,  not  the  image  itself.  If  we  are  only  trying  to 
measure  binary  star  separations  we  may  have  all  the  information  we  need,  but  high-resolution  imaging 
does  not  fall  out  immediately.  We  can  take  V  I  F(m,v)  |  ^  and  obtain  the  correct  magnitude  of  F  (u,v) 
but  we  are  lacking  the  phases. 

The  phase  retrieval  problem  will  be  discussed  in  more  detail  below,  but  suffice  it  to  say  that  there 
exist  practical  methods  for  reconstructing  the  phase  of  F  (u,v)  in  many  cases.  This  is  an  area  of  current 
research,  but  major  steps  have  been  taken  in  the  last  ten  years.  Coupling  Labeyrie’s  method  for  estimat¬ 
ing  I  F (u,v)  I  with  one  or  more  of  the  phase  reconstruction  algorithms  yields  a  complete  high- 
resolution  imaging  capability. 
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3.  STATISTICAL  OPTICS 


3.1  INTRODUCTION 

The  optical  theory  required  by  this  report  is  well  presented  in  the  book  by  Goodman.^  We  will 
extract  some  information  from  it  in  roughly  the  same  order  Goodman  presents  it.  Goodman  devotes  two 
chapters  to  summarizing  the  theory  of  random  processes.  A  more  expanded  treatment  of  this  topic  may 
be  found  in  Papoulis^  and  similar  texts.  We  will  assume  that  the  reader  is  familiar  with  this  material. 

An  additional  topic  of  interest  is  a  short  review  of  the  theory  of  analytic  signals.  Appendix  A  will 
discuss  this  subject  in  as  much  detail  as  required  here. 

The  first  topic  of  the  next  subsection  will  be  first-order  properties  of  light,  i.  e.,  light  properties  at  a 
single  instant  in  time.  We  will  discuss  propagation  of  light  waves  of  various  bandwidths  using  the 
Huygens-Fresnel  principle.  We  will  define  the  instantaneous  intensity  (that  quantity  to  which  optical 
detectors  respond). 

At  this  point,  we  will  consider  second-order  properties  of  light  in  the  form  of  coherence  functions. 
These  include  temporal  and  spatial  coherence  functions,  the  mutual  coherence  function,  and  mutual 
intensity.  This  material  is  needed  to  understand  the  operation  of  interferometers,  the  sensing  elements  of 
adaptive  optics. 

Propagation  of  mutual  intensity  from  an  incoherent  source  is  described  by  the  Van  Cittert-Zemike 
theorem,  the  basis  for  speckle  imaging.  This  theorem  states  that  the  mutual  intensity  is  the  two- 
dimensional  Fourier  transform  of  the  intensity  distribution  across  the  source.  This  is  reminiscent  of  the 
relationship  of  the  illumination  function  of  an  antenna  and  the  far  field,  but  the  quantities  are  entirely 
different. 

Other  material  falling  in  the  category  of  statistical  optics  relates  directly  to  succeeding  major  sec¬ 
tions  of  this  report,  and  will  be  covered  in  those  sections. 

3.2  PROPAGATION  OF  LIGHT 

3.2.1  Propagation  of  Monochromatic  Light 

Consider  the  scalar  amplitude  of  one  polarization  component  of  the  electrical  or  magnetic  field  of 
monochromatic  light.  If  P  represents  a  position  in  space  and  t  a  point  in  time,  we  will  denote  this  ampli¬ 
tude  u  The  analytic  signal  representation  of  u{P,()  is 

u(P,r)  =  U(P,v)e-^'2’^'  (3.1) 

where  v  is  the  frequency  of  the  wave  and  U(P,  v)  is  its  amplitude.  Let  such  a  wave  be  incident  from  the 
left  on  the  infinite  surface  Z  shown  in  Figure  3.1. 

We  are  interested  in  the  phasor  amplitude  of  the  field  at  the  point  Pq  in  terms  of  the  field  on  L.  We 
will  use  the  approximation  known  as  the  Huygens-Fresnel  principle.  If  c  is  the  sp)eed  of  light,  X  =  -^  is 
the  wavelength,  x(0)  is  an  obliquity  factor  with  x(0)  =  1  and  0  <  x(0)  ^  1,  and  r^X, 
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Figure  3.L  Propagation  geometry. 
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Note  from  Figure  3.1  that  0  is  the  angle  between  the  line  connecting  Pq  and  P  i  and  the  normal  to  Z.  In 
words,  each  point  on  Z  acts  as  a  secondary  source  of  spherical  waves  radiating  with  directional  ampli¬ 
tude  pattern  x(0).  More  detail  on  the  Huygens-Fresnel  principle  may  be  obtained  from  any  text  on 
optics. 

3.2.2  Propagation  of  Nonmonochromatic  Light 

Since  monochromatic  light  is  of  limited  interest  to  us,  we  want  to  generalize  our  expression  for 
light  propagation.  Let  u{P,t)  be  a  wave  whose  analytic  representation  is  denoted  u(F,r).  It  is  possible  to 
apply  the  Huygens-Fresnel  principle  inside  of  an  integral  over  frequency  to  obtain  an  expression  for  gen¬ 
eral  propagation.  The  main  complication  is  that  since  u{P,t)  may  not  be  Fourier  transformable,  a  trun¬ 
cated  interval  is  used  in  the  derivation  and  a  limit  is  later  taken.  After  a  little  manipulation,  we  obtain 


—  iiCP.  t  -rlr\ 


(3.3) 


3.2  J  Propagation  of  Narrowband  Light 

If  the  light  bandwidth  Av  is  much  smaller  than  the  center  frequency  v,  the  propagation  condition 
may  be  simplified  to 


(3.4) 


6 


3.3  OPTICAL  COHERENCE 


3.3.1  Temporal  Coherence 

If  u{Pj)  is  the  analytic  signal  corresponding  to  an  optical  wave,  consider  the  relationship  between 
u(P,r)  and  u(P,r+T).  u{Pj)  has  some  bandwidth  Av,  and  hence  the  complex  envelope  I  u{Pj)  I  is 
varying  at  a  rate  determined  by  Av.  If  x  1/Av  we  expect  u(P,r)  and  u(p,r+x)  to  be  highly  correlated. 

At  this  point,  we  will  digress  with  a  word  about  detectors.  Optical  detectors  respond  to  the  inten¬ 
sity  of  the  optical  wave,  e.  g.,  the  square  of  the  magnitude  of  the  analytic  signal  representation.  We 
know  the  magnitude  of  the  analytic  signal  as  the  complex  envelope.  The  only  thing  to  keep  in  mind  is 
that  most  detectors  average  over  a  considerable  amount  of  time.  We  are  really  recording  some  average 
of  the  square  of  the  complex  envelope. 

TTie  canonical  experiment  in  temporal  coherence  was  p>erformed  by  Michelson  in  the  19th  century. 
In  the  Michelson  interferometer,  a  point  source  of  light  is  split  by  a  beam  splitter  and  reflected  from  two 
different  perpendicular  mirrors,  producing  a  selected  path  length  difference.  A  detector  records  the 
resulting  average  intensity. 

Mathematically,  if  we  let  the  path  difference  be  /?,  the  time  delay  is  x  =  /?/r.  The  measured  inten¬ 
sity  will  be 

/(t)  =  <lu(0  +  u(r+t)|2>  (3.5) 

where  the  angle  brackets  denote  a  time  average.  This  will  look  something  like  Figure  3.2.  For  small 
path  differences  (time  delays)  everything  is  in  phase  and  reinforces.  As  the  delay  increases,  we  see  peri¬ 
odic  behavior  at  the  average  wavelength  X,  but  for  larger  delays,  the  multiple  frequencies  present  in  the 
wave  cause  dispersion  and  eventually  result  in  a  constant  value  with  no  fringes.  There  are  fringes 
present  at  each  individual  frequency,  but  they  add  up  in  random  phase  to  a  constant  value. 

The  exact  behavior  of  a  given  optical  wave  depends  on  the  bandwidth  of  that  wave,  or,  more  gener¬ 
ally,  on  its  spectral  density.  Much  more  can  be  said  about  this  subject  but  we  will  not  need  it  for  our 
development. 


If  we  expand  Equation  (3.5)  and  define  the  average  intensity 

Io={  \  U(r)  1 

(3.6) 

and  the  self-coherence  function 

r(T)s(u(r-i-T)  u*(0> 

(3.7) 

we  obtain 

/(T)  =  2|/o-HRe  [Rt)]  ). 

(3.8) 
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Figure  3.2.  Typical  fringe  pattern  with  envelope. 

3.3.2  Spatial  Coherence 

33.2.1  Young’s  Experiment 

Every  real  light  source  has  some  bandwidth,  so  the  temporal  coherence  effects  discussed  above  are 
important.  In  addition,  every  real  source  has  some  finite  size,  so  the  point  source  approximation  is  never 
completely  valid.  When  we  consider  extended  sources,  we  need  to  consider  both  spatial  and  temporal 
coherence. 

In  our  discussion  of  temporal  coherence,  we  considered  the  analytic  signals  u(P,t)  and  u(P,r+T). 
For  X  smaller  than  the  correlation  time,  the  signals  were  correlated.  In  the  case  of  spatial  coherence  we 
want  to  consider  the  analytic  signals  u(P  j ,/)  and  u(P2>0-  We  expect  that  for  point  separations  smaller 
than  the  correlation  distance,  the  signals  will  be  correlated. 

It  is  possible  to  work  with  close  approximations  to  point  sources  and  get  results  in  good  agreement 
with  the  theory  of  the  previous  subsection.  Such  light  has  a  spatial  correlation  function  which  is  an 
impulse,  e.  g.,  its  spatial  frequency  spectrum  is  white. 

Unfortunately,  it  is  not  practical  to  obtain  light  which  has  an  impulse  for  its  time  correlation  func¬ 
tion  (white  spectrum).  Hence  we  are  unable  to  study  the  effects  of  spatial  correlation  in  isolation;  we 
must  consider  spatial  and  temporal  correlation  together. 

The  canonical  experiment  in  spatial  and  temporal  correlation  was  performed  by  Thomas  Young. 
As  shown  in  Figure  3.3,  light  from  an  extended  source  passes  through  two  pinholes  and  then  forms  an 
interference  pattern  on  a  viewing  screen. 
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Figure  3.3.  Young* s  intcfference  experiment. 


If  we  assume  that  the  lime  delay  difference  {r2-r\)lc  is  much  less  than  the  coherence  time  of  the 
light  source,  interference  fringes  will  appear  on  the  viewing  screen.  The  depth  of  modulation  will 
depend  on  the  degree  of  correlation  between  the  light  waves  from  the  two  points.  We  express  the  inten¬ 
sity  at  point  Q  as 

/(0)  =  <u*(0,Ou(0,O>.  (3.9) 


3  J.2.2  Narrowband  Case 

If  the  light  is  narrowband,  and  the  pinholes  are  not  too  large,  we  can  write  the  analytic  signal  at  Q 
in  terms  of  the  signals  at  P  \  and 

\x{QJ)^KMPy,t-r^lc)-^K2^^{P2^t^r2lc).  (3.10) 

We  can  employ  Equation  (3.4)  to  perform  the  unusual  feat  of  integrating  over  pinholes  and  find  that 

f  f  ^(0/) 

A:/=  /  =  1,2.  (3.11) 

p,  j'^i 

It  is  not  our  purpose  to  impart  great  insight  into  the  results  of  Young’s  experiment  and  all  the  pos¬ 
sible  fringe  patterns  that  result,  but  we  do  want  to  see  how  even  in  such  a  simple  experiment  certain 
multi-dimensional  quantities  appear.  It  will  be  important  for  our  later  analysis  to  understand  the  basic 
definitions  and  units  of  some  of  these  quantities. 

To  that  end,  define 

•  r‘\Q)=  \Ki\^{\u(Pi,t-r^/c)\^),  i  =  \,2  (3.12) 

the  individual  intensities  and 

r,2(x)  =  <u(P, ,r  X)  u‘(P2,0>  (3.13) 
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the  mutual  coherence  function.  Substituting  these  expressions  into  Equations  (3.10)  and  (3.9),  and  not¬ 
ing  that  the  Kj  are  purely  imaginary,  we  obtain 

I(Q)  =  f^HQ)  +  /<‘)((2)  +  2K  Re[r,2(-^^^)]  (3.14) 

c 

where 

K  =  K^K\^K\K2.  (3.15) 

We  are  studying  the  intensity  of  an  optical  wave,  which  has  units  of  power.  We  have  decomposed 
the  intensity  resulting  from  Young’s  experiment  into  three  parts:  the  intensities  contributed  by  the  dif¬ 
fraction  patterns  of  the  first  and  second  pinholes  and  the  intensity  resulting  from  mutual  coherence.  The 
actual  pattern  of  intensity  is  rather  complicated  and  reflects  both  temporal  and  spatial  coherence. 

Since  the  pinholes  are  small,  the  diffraction  patterns  are  essentially  constant.  Hence  they  contrib¬ 
ute  a  constant  background  upon  which  is  superimposed  a  fringe  pattern  with  slowly  varying  envelope 
and  phase.  Where  the  path  length  difference  is  zero  or  small,  the  fringe  pattern  is  primarily  attributable 
to  temporal  coherence  effects.  When  the  path  length  difference  becomes  large,  the  shrinking  of  the 
fringe  envelope  is  primarily  attributable  to  spatial  coherence  effects. 


3.3.2.3  Quasi-monochromatic  Case 


We  want  to  consider  an  even  more  special  case  in  this  subsection.  In  addition  to  the  narrowband 
assumption  in  which  Av  •«  v  we  will  here  assume  that  the  coherence  length  of  the  light  is  much 
greater  than  the  maximum  path  length  difference.  Referring  to  Figure  3.3, 

(,-i  +  ri/)-(r2 -H  r2') 

- (3.16) 

c 

With  this  assumption,  the  fringe  contrast  is  constant  over  the  observation  region  of  interest  and  the 
mutual  coherence  function  can  be  written  as 


r,2(T)  =  J,2^-''’^^  :  (3.17) 

where  J12  =  ri2(0)  is  the  mutual  intensity  of  the  light  at  pinholes  P\  and  P2*  ^  this  formulation,  J 12  is 
the  phasor  amplitude  of  a  constant  spatial  sinusoidal  fringe  across  the  observation  region. 


33.3  Propagation  of  Mutual  Coherence 

As  an  optical  wave  propagates  through  space,  its  detailed  structure  changes.  Similarly,  its  mutual 
coherence  function  changes,  so  we  can  say  that  mutual  coherence  “propagates.”  Things  get  a  little  more 
complicated  when  we  consider  the  mathematical  representation  of  propagation  of  mutual  coherence. 

As  can  be  seen  with  reference  to  Figure  3.4,  the  dimensionality  of  the  functions  increases.  Instead 
of  two  pinholes  P  j  and  P2»  have  two  sets  of  points  Pj  and  Qf.  We  assume  that  we  know  the  mutual 
coherence  function  r(P  1  ,^2^^)  on  the  surface  Zi ,  that  is,  we  know  the  results  of  Young’s  interference 
experiment  for  all  possible  pairs  of  pinholes  P,  .  We  wish  to  use  this  information  to  determine  the  mutual 
coherence  function  on  the  surface  2^,  that  is,  predict  the  results  of  Young’s  interference 

experiment  for  all  possible  pairs  of  pinholes  (2/. 
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Figure  3.4.  Mutual  coherence  propagation  geometry. 


Modifying  our  notation  slightly,  in  Equation  (3.13)  we  defined  the  mutual  coherence  function  on 

1.2  as 

r(!2 1  >02;'^)  =  (u((2 1  +'r)  u  ((22>o)-  (3.18) 


We  can  use  Equation  (3.4)  to  express  the  fields  on  I2  in  terms  of  the  fields  on  Lj  in  the  case  of  nar¬ 
rowband  light: 


f  f u(Pi,r  +  x-rxlc) 

u(0,,r+T)=  JJ - ; - X(0i)^5, 

li  j^rx 

,  e  e-\l  {P 2,t -rile) 

U  (02,0=  JJ - ^ - —X(.^2)dS2 


jh-2 


(3.19) 


Substituting  (3.19)  into  (3.18),  with  a  little  manipulation  and  another  application  of  (3.19)  we 
obtain  the  moderately  formidable 


r(A,,02;O=  JJ  JJr(P,,P2;T  +  -^ 


X(9i)  X(92) 
h-x  ^2 


dSx  as  2 . 


(3.20) 


We  will  not  need  to  spend  any  time  contemplating  the  inner  nature  of  quadruple  integrals;  we  will  imme¬ 
diately  specialize  the  result  and  get  things  to  more  manageable  proportions.  Our  goal  is  to  provide  back¬ 
ground  helpful  in  understanding  the  important  Van  Cittert-Zemike  theorem,  and  we  are  almost  there. 


33.4  Propagation  of  Mutual  Intensity 

If  the  quasi-monochromatic  condition  of  Equation  (3.16)  is  satisfied  the  mutual  intensity 

J(0 1, 02)  =  n0 1,02:0) 

is  of  interest.  Applying  Equation  (3.17)  to  Equation  (3.20),  we  obtain 
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(3.21) 


-/27C{r')-r  1 ) 

-  -  x"  -  xcQi)  x(e2) 


dSi  dS2, 


(3.22) 


MQi,Q2)=  11  11  J(^l,^2)^ 

ii  i|  X/'i  Xr2 

the  propagation  law  for  mutual  intensity.  The  intensity  distribution  on  the  surface  is  obtained  by 
letting  Q\  Q2  in  Equation  (3.22). 


3.4  THE  VAN  CITTERT-ZERNIKE  THEOREM 
3.4.1  Coherent  and  Incoherent  Sources 


We  will  next  briefly  discuss  coherent  and  incoherent  wavefields.  It  seems  logical  to  call  optical 
waveforms  at  points  P 1  and  ^2  subject  to  time  delay  T  fully  coherent  if  the  magnitude  of  the  normalized 
coherence  function  is  unity 

I  Yi2(t)  I  =  1  (3.23) 


where 

ri2(T) 

Vrn(0)r22(0)  • 


(3.24) 


Y  ranges  from  0  to  1  in  the  same  way  as  a  correlation  coefficient.  Since  only  monochromatic  waves  sat¬ 
isfy  (3.23),  it  has  been  found  useful  to  call  a  wavefield  fully  coherent  if,  for  every  pair  of  points  P 1  and 
P2,  there  exists  some  delay  T  such  that  Equation  (3.24)  is  satisfied. 

Similarly,  it  seems  logical  to  define  an  incoherent  field  by  requiring 

|ri2(T)l=0  Pi^P2  (3.25) 


but  this  definition  is  even  less  useful  than  the  first  attempt  at  defining  fully  coherent  wavefronts  since  it 
does  not  allow  any  propagation  at  all!  Essentially  it  describes  a  nonradiating  source.  The  closest  we  can 
come  to  incoherence  with  a  radiating  wave  turns  out  to  be  a  fairly  messy  expression  involving  Bessel 
functions.  However,  in  practice  as  long  as  the  optical  system  resolution  is  much  less  than  X,  we  can 
approximate  the  mutual  intensity  of  an  incoherent  source  by 

J{P  I, P  2)  =  — HP  Ay)  (3.26) 

7t 

where  5(.,.)  is  a  two-dimensional  Dirac  delta  function  and  Ax  and  Ay  are  the  differences  in  the  x  and  y 
coordinates  of  points  P 1  and  /’2- 


3.4.2  Outline  of  Mathematical  Derivation 

The  Van  Cittert-Zemike  theorem  describes  the  mutual  intensity  function  produced  by  an  incoher¬ 
ent  source.  Such  sources  are  of  primary  concern  to  us  in  the  passive  surveillance  situation,  as  discussed 
previously.  If  we  start  with  Equation  (3.22),  the  propagation  law,  and  substitute  in  Equation  (3.26)  for 
the  mutual  intensity,  we  nicely  simplify  things  to 
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(3.27) 


-j2n{r2-rO 

J(01,02)  =  -  JJ  nP^e  ^ 

"  I 


X(9i)  X(92) 
7ir\  7j‘2 


dS 


where  I  is  the  source. 

This  result  is  useful  for  quasi-monochromaiic  incoherent  sources  only,  but  we  need  to  sf>ecialize  it 
still  further.  Assume  that  the  distance  r  between  the  source  and  the  observation  region  is  much  greater 
than  the  size  of  either  of  them.  Then  r]r2  -2^,  and  since  the  angles  are  small,  X(0/)^  ^hen 

obtain  the  simplification 

-j2n{r2-ri) 

J((2i,02)  =  -T  n  ^  (3.28) 

Jt-  I 


To  define  the  coordinate  systems,  let  the  source  and  observation  planes  be  separated  in  the  z  dimension, 
let  the  perpendicular  coordinates  in  the  observation  plane  be  {x,y),  and  let  the  coordinates  in  the  source 
plane  be  (^,1)).  Substituting  these  coordinates  for  the  r,  and  using  the  approximations,  with  a  little 
manipulation  we  arrive  at  the  expression 


J(.v,,y,;.V2,V2) 


p-jv  ff  p-j2n(^xi,  + &yr\) 

^JJ  I(^,T1)^^ - - - dS 

nz^  I  X 


(3.29) 


where  the  phase  factor  Vj/  is  given  by 


+  >’2  )  -  (^1 


(3.30) 


In  words,  the  mutual  intensity  function  at  the  observation  plane  is  the  two-dimensional  Fourier 
transform  of  the  source  intensity  multiplied  by  a  phase  factor.  There  are  conditions  under  which  the 
phase  factor  is  approximately  unity  and  may  be  omitted.  In  our  situation,  we  will  not  need  to  worry 
about  this  factor. 

Note  that  mutual  intensity  is  more  general  than  image  intensity;  when  discussing  the  latter, 
J(.v ]  ,y  I  j^  2^y2)  becomes  simply  JCv,y)  as  the  two  points  converge  to  one. 


3.5  PROPAGATION  THROUGH  TRANSMITTING  OBJECTS 

Goodman"^  presents  the  general  theory  which  enables  one  to  compute  image  intensity  when  par¬ 
tially  coherent  light  is  transmitted  through  objects.  It  is  then  easy  to  consider  the  special  cases  of  total 
coherence  and  incoherence.  Since  things  get  rather  complicated  and  we  do  not  need  to  consider  coherent 
sources,  we  will  just  present  the  results  we  need  without  summarizing  any  of  the  proofs. 

We  will  first  describe  the  intensity  Fourier  transform  property  of  lenses.  This  will  provide  the 
mathematical  framework  for  understanding  imaging  with  incoherent  light.  We  will  then  define  the  opti¬ 
cal  transfer  function.  This  will  be  important  in  characterizing  the  atmosphere. 
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3.5.1  The  Effect  of  a  Lens  on  Image  Intensity 

It  is  well  known  that  a  lens  acts  as  a  Fourier  transformer  for  coherent  light;  the  complex  fields  in 
the  focal  planes  are  a  Fourier  transform  pair.  We,  however,  are  concerned  with  incoherent  light,  which 
leads  us  to  consider  the  propagation  of  mutual  coherence  through  a  lens. 

If  we  restrict  ourselves  to  quasi-monochromatic  light,  it  can  be  shown  that  the  mutual  intensity 
functions  in  the  focal  planes  of  a  thin  positive  lens  are  a  four-dimensional  Fourier  transform  pair.  We  are 
most  interested  in  the  intensity  in  the  back  focal  plane  which  results  from  an  incoherent  source  in  the 
front  focal  plane.  In  this  case,  both  functions  reduce  from  four  dimensions  to  two  and  we  have  a  two- 
dimensional  Fourier  transform  relationship. 

We  will  not  be  studying  this  derivation  or  performing  any  calculations  with  it,  so  we  will  not  show 
the  relevant  equations.  The  theorem  we  wish  to  remember  is  simply  that  when  we  have  quasi- 
monochromatic  incoherent  light  at  the  front  focal  plane  of  a  thin  lens,  the  intensity  distribution  in  the 
back  focal  plane  is  the  two-dimensional  Fourier  transform  of  the  source  mutual  intensity. 

This  theorem  complements  the  Van  Cittert-Z^mike  theorem  discussed  above.  Assume  an  astro¬ 
nomical  object  at  a  great  distance  is  illuminated  by  (or  emits)  incoherent  light  which  is  either  inherently 
quasi-monochromatic  or  is  later  filtered  to  satisfy  this  condition.  By  Van  Cittert-Z^mike  the  mutual 
intensity  is  the  Fourier  transform  of  the  source  intensity.  If  the  light  then  enters  a  suitable  lens,  the  back 
focal  plane  intensity  will  be  the  Fourier  transform  of  the  mutual  intensity.  Since  the  Fourier  transform 
and  the  inverse  Fourier  transform  differ  only  in  scale  factor  and  phase,  we  see  an  image  of  the  source, 
which  probably  won’t  surprise  anyone  who  believes  in  telescopes. 

One  missing  factor  in  this  scenario  is  the  atmosphere,  which  degrades  the  image  so  that  it  is  not  a 
replica  of  the  source  intensity  distribution.  In  addition,  the  finite  size  of  the  lens  means  that  we  are 
selecting  a  finite  number  of  spatial  frequencies,  and  hence  there  is  a  high-frequency  cutoff  determined  by 
the  size  of  the  lens  pupil. 

3.5.2  The  Optical  Transfer  Function 

There  are  certain  conditions  under  which  a  transmitting  object  can  be  characterized  by  an  optical 
transfer  function  (OTF).  Such  an  object  could  be  a  lens  or  other  imaging  system,  the  atmosphere,  or 
some  combination. 

First,  the  illumination  must  be  incoherent  in  the  sense  of  Equation  (3.26).  We  will  find  that  an 
incoherent  system  is  linear  in  intensity,  while  a  coherent  system  is  linear  in  complex  amplitude. 

The  second  requirement  for  the  OTF  to  exist  is  that  our  system  must  be  isoplanatic.  The  spatial 
analog  of  the  stationarity  property  of  time  processes  goes  by  a  number  of  names,  including  shift- 
invariance,  space-invariance,  and  isoplanaticity.  By  any  name,  it  loosely  means  that  the  statistical  pro¬ 
perties  of  a  random  process  do  not  depend  on  the  spatial  coordinates  at  which  they  are  measured.  In 
practice  this  is  true  if  the  object  under  examination  does  not  subtend  a  very  large  angle. 

Under  these  assumptions,  and  assuming  quasi-monochromaticity,  if  we  denote  the  object  intensity 
distribution  o(x,y)  and  the  transmitting  object  point  spread  function  (really,  intensity  spread  function) 
t{x,y),  we  find  that  the  image  intensity  is 

/  (j:,y)  =  o  {x,y)  *  t  (A-,y).  (3.31) 
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This  convolution  relationship  is  a  product  in  the  frequency  domain 

/(w,v)  =  0(w,v)7'(w,v).  (3.32) 

Note  that  we  have  here  departed  from  Goodmein’s  notation,  emticipating  future  discussions  of  speckle 
image  processing. 

Intuitively,  when  a  distemt  object  is  imaged  through  the  atmosphere,  we  apply  first  Van  Cittert- 
Zemike  and  then  the  thin  lens  relationship.  The  atmosphere  OTF  (9(w,v)  is  acting  on  the  mutual  inten- 
sity  function,  the  Fourier  tremsform  of  the  source  intensity.  If  we  insist  on  considering  things  in  the 
intensity  domain,  a  convolution  with  the  point  spread  function  o  {x,y)  is  necessar>'. 

The  OTF  is  a  complex  quantity;  the  magnitude  of  the  OTF  is  known  as  the  modulation  transfer 
function  (MTF)  and  will  be  importemt  in  our  later  discussions. 

3.6  THE  INTERFEROMETRIC  VIEW  OF  IMAGING 

Recall  Young’s  experiment,  in  which  light  passes  through  two  pinholes  and  forms  an  interference 
pattern  on  a  viewing  screen  as  in  Figure  3.3.  If  the  pinholes  are  moved  close  together,  a  one-wavelength 
path  difference  spans  more  distance  on  the  screen,  so  the  period  of  the  interference  pattern  is  large. 
Greater  pinhole  spacing  produces  interference  patterns  with  smaller  period  or  higher  spatial  frequency. 
The  pinholes  are  taken  to  be  small  compared  to  the  coherence  distemce  of  the  light  source. 

The  exit  pupil  of  an  imaging  system  can  be  assumed  to  be  composed  of  many  pinholes.^’ Each 
possible  pair  of  these  imaginary'  pinholes  generates  an  interference  pattern,  and  we  can  consider  the 
whole  image  to  be  built  up  of  a  number  of  sinusoidal  fringes.  The  spatial  frequencies  and  are 
related  to  the  pinhole  separations  Ax  and  Ay  by 

Av  = 

Ay  =  ^v,..  (3.33) 

If  the  object  is  at  a  great  distance  z  it  is  convenient  to  divide  both  sides  of  these  equations  by  z,  express¬ 
ing  the  separation  in  radians  rather  than  meters  and  the  spatial  frequencies  in  cycles/rad  rather  than 
cycles/m. 

The  amplitude  and  phase  of  the  fringe  pattern  from  a  given  pair  of  pinholes  located  at  (Jti,y  j)  and 
(A'2,y2)  determined  by  the  amplitude  and  phase  of  the  mutual  intensity  leaving  the  exit  pupil 
Since  there  are  many  pinhole  pairs  with  any  given  spacing,  the  resultant  fringe  pattern 
at  a  given  spatial  frequency  is  calculated  by  summing  over  all  such  pairs.  If  the  optics  are  perfect,  cind 
assuming  incoherent  illumination  cind  Van  Cittert-Zemike,  the  mutual  intensity  is  of  the  form  Jp{Ax,Ay) 
and  all  pinhole  pairs  with  a  given  spacing  add  coherently  to  form  the  image. 

Consider  the  transfer  function  of  an  ideal  circular  lens.  It  is  easy  to  see  that  the  higher  spatial  fre¬ 
quencies,  corresponding  to  greater  point  spacings,  are  attenuated  relative  to  the  lower  frequencies,  since 
there  are  fewer  point  pairs  in  the  aperture  with  large  spacings.  At  left  in  Figure  3.5  are  representative 
point  pairs  with  identical  spacing.  The  number  of  such  pairs  relative  to  the  number  of  pairs  with  zero 
spacing  is  given  by  the  shaded  area  in  the  right  diagram  divided  by  the  total  aperture  size.  This  is  mathe¬ 
matically  an  autocorrelation  of  the  circular  lens  pupil. 
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(a)  Typical  Point  Pairs 


(b)  Weighting  Factor 


Figure  3.5.  Interferometric  imaging  model. 
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4.  ATMOSPHERIC  TURBULENCE 


4.1  INTRODUCTION 

The  general  theor>'  of  light  propagation  through  atmospheric  turbulence  was  introduced  in  the 
important  book  by  Tatarskii  This  translation  from  the  Russian  drew  on  a  lot  of  Russian  work  and  has 
been  the  basis  for  most  of  the  work  in  this  area  since.  The  mathematics  is  not  simple,  so  we  will  again 
adopt  the  approach  of  outlining  the  important  results  in  an  effort  to  gain  intuitive  understanding. 

After  Tatarski,  we  wish  to  discuss  briefly  the  papers  of  Hufnagel  and  Stanley^  and  Fried.^  Another 
major  step  in  this  area  was  made  by  Korff*^  but  we  will  consider  it  in  our  discussion  of  sp)eckle  imaging 
since  it  arose  in  that  context.  All  these  authors  rely  on  the  approach  of  Tatarski,  and  use  similar  notation. 
We  will  describe  the  assumptions  made  by  Tatarski  and  his  successors  about  the  atmosphere,  as  well  as 
the  effect  of  those  models  on  the  amplitude  and  phase  of  transmitted  light.  We  will  later  apply  these 
results  to  the  analysis  of  adaptive  optics  and  speckle  imaging. 

Not  all  of  Tatarski ’s  assumptions  seem  transparent  on  first  reading;  in  particular  it  is  sometimes 
hard  to  determine  typical  ranges  of  certain  parameters.  Goodman‘S  puts  Tatarski’s  work  in  slightly  more 
modem  context  and  is  recommended  as  an  introduction. 

4.2  STATISTICAL  REPRESENTATION  OF  ATMOSPHERIC  PARAMETERS 
4.2.1  Structure  Functions 

It  is  convenient  when  dealing  with  stationary  random  functions  to  use  the  mean,  the  autocorrela¬ 
tion  function  R  (i),  and  the  spectral  density  S  (/),  the  latter  two  being  of  course  a  Fourier  transform  pair. 
When  a  process  is  nonstationary',  we  must  work  with  the  more  general  mean,  a  function  of  time,  and  the 
more  general  autocorrelation  function  R  {r,s). 

In  meteorology,  we  often  deal  with  quantities  which  are  nonstationary  such  as  air  pressure  and 
temperature.  About  50  years  ago  Kolmogorov  introduced  a  means  of  mathematically  separating  very 
slow  variation  in  the  mean  from  the  more  rapid  variation  for  any  random  process.  Instead  of  using  the 
absolute  value  of  some  parameter  /  (/)  he  considered  the  incremental  changes 

F,(/)^/(/+t)  -/(/).  (4.1) 

It  is  intuitively  obvious  that  if  T  is  chosen  appropriately,  we  may  end  up  with  an  essentially  stationary 
process  which  reflects  all  fluctuations  whose  period  is  less  than  x. 

The  next  obvious  step  is  to  consider  the  autocorrelation  function  of  the  process  Fj.  It  turns  out, 
however,  as  discussed  in  Appendix  B,  that  this  autocorrelation  function  can  be  expressed  in  terms  of  so- 
called  structure  function  If  the  structure  function  may  be  represented  in  the  form 

Dfix)  =  |/(r  +  T)-/(0f  (4.2) 

we  say  that  the  process  has  stationary  increments.  This  is  the  model  of  interest  to  us. 
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4.2.2  Random  Fields 


The  multi-dimensional  analog  of  the  random  process  is  called  a  random  field.  Notationally,  fol¬ 
lowing  Tatarski,  in  the  three-dimensional  case  define  the  vector =  (^,^,"2;).  The  spatial  autocorrela¬ 
tion  function  of  a  random  field  takes  the  general  form  R  j  ,7^2)-  l^e  random  field  is  homogeneous  and 
isotropic,  this  may  be  written  as  /?(  IT^i  -1^2  I  )•  In  this  latter  case,  things  reduce  to  a  single  dimension; 
in  particular,  the  spectral  density  exists  as  a  function  of  one  variable. 

It  may  seem  like  quite  a  simplification  to  model  the  atmosphere  as  homogeneous  and  isotropic.  In 
fact,  the  large  scale  effects  which  cause  slow  variation  in  means  of  atmospheric  parameters  also  result  in 
inhomogeneities  and  anisotropy.  The  solution  is  the  same;  the  use  of  structure  functions  enables  us  to 
obtain  approximately  homogeneous,  isotropic,  stationary  random  fields. 

4.3  PHYSICS  OF  TURBULENCE 

Turbulence  exists  on  a  wide  range  of  scales.  The  smallest  of  these  is  in  the  range  of  several  milli¬ 
meters,  and  the  largest  is  on  the  order  of  tens  of  meters.  Kolmogorov  denoted  these  the  “inner”  and 
“outer”  scales  of  the  turbulence. 

The  atmospheric  quantity  of  most  interest  to  us  is  the  refractive  index.  This  depends  on  tempera¬ 
ture,  pressure,  and  wavelength,  but  the  fluctuations  depend  mostly  on  temperature.  Define  the 
wavenumber  vector (Kj,Ky,K-),  a  vector  frequency  measured  in  radians  per  meter.  In  the  isotropic 
case  we  can  work  with  K  =  •  Fof  ^  consistent  with  the  two  turbulence  scales  mentioned 

above  Kolmogorov  gave  the  form  of  the  power  spectral  density  of  refractive  fluctuations 

<D„(k)  =  0.033  k-"^3  (4.3) 

In  this  expression,  is  called  the  structure  constant  of  the  refractive  index  fluctuations.  When  we 
couch  things  in  terms  of  structure  functions, 

D„{r)  =  Clr^'\  (4.4) 

The  numerical  value  of  ranges  from  10“'^  to  10“’^ 

Kolmogorov’s  work  is  rather  old,  but  like  politicians  it  is  easier  to  criticize  than  to  replace.  One 
new  area  of  study  which  is  relevant  to  turbulence  is  the  area  of  nonlinear  dynamics,  often  called  chaos 
theory.”  Kolmogorov  was  aware  that  the  above  expression  is  not  applicable  when  energy  dissipation 
becomes  too  significant  (for  small  scale  sizes)  but  more  is  understood  about  the  phenomenon  today. 
Clearly  there  is  merit  to  Kolmogorov’s  work  since  it  is  supported  by  experiment  to  useful  extent,  but 
future  advances  may  give  us  more  accurate  analyses. 

4.4  OPTICAL  PROPAGATION  THROUGH  TURBULENCE 
4.4.1  Amplitude  and  Phase  Fluctuations 

It  is  beyond  the  scope  of  this  report  to  reproduce  much  of  Tatarski ’s  analysis  of  optical  propagation 
through  turbulence.  We  will  not  consider  Maxwell’s  equations  or  a  host  of  other  relationships  and 
approximations  which  are  relevant.  Our  goal  is  just  to  present  the  end  results,  with  as  much  physical 
meaning  as  possible. 


18 


Tatarski  moves  through  a  series  of  successive  analyses  of  progressively  more  realistic  scenarios. 
Although  it  might  build  character  for  us  to  recapitulate  some  of  this,  we  will  instead  skip  to  the  most  rel¬ 
evant  case.  Here  we  assume  that  the  structure  constant  Cl  varies  slowly,  with  significant  changes  only 
occurring  for  path  lengths  on  the  order  of  the  outer  scale  of  turbulence. 

We  will  present  expressions  for  the  log-amplitude  and  phase  of  the  fluctuation.  There  are  a  couple 
of  reasons  why  we  consider  the  log-amplitude  rather  than  just  the  amplitude,  but  suffice  it  to  say  here 
that  for  weak  fluctuations  (good  seeing)  experiment  suggests  that  amplitude  fluctuations  are  log- 
normally  distributed.  Equivalently,  the  log-amplitude  is  normally  distributed.  Hence  the  log-amplitude 
is  the  natural  quantity  of  interest. 

Notationally,  let  the  refractive  index  be  represented  as 

n(i^  =  no(i^  +  n^lj^  (4.5) 

where  the  constant  term  nQ(i^-  \  and  the  fluctuating  term  1.  Let  the  phase  of  the  wave  be 

5  =  5o  +  5]  and  the  amplitude  A  satisfy  logi4  =logy4o  +  X-  Tatarski  gives  an  expression  for  the 
log-amplitude  fluctuations 

_  L 

(4.6) 

0 

We  won’t  worry  about  the  constant  since  we  won’t  be  doing  any  numerical  calculations  with  this 
expression. 

The  phase  of  the  fluctuation  is  generally  assumed  normally  distributed.  Rather  than  considering 
the  fluctuations  directly,  Tatarski  derives  an  expression  for  the  structure  function  of  the  phase 

L 

D5(p)  =  /f5p5'3  (4.7) 

0 

Again,  the  exact  value  of  the  constant  is  not  of  interest. 

To  interpret  Equations  (4.6)  and  (4.7),  consider  the  following  facts.  First,  amplitude  fluctuations 
are  much  less  significant  than  phase  fluctuations;  indeed  adaptive  optics  and  speckle  imaging  would  not 
work  in  their  present  forms  if  this  was  not  the  case.  Second,  since  the  density  of  air  decreases  so  rapidly 
with  elevation,  refractive  index  changes  tend  to  be  greatest  near  the  ground. 

Looking  at  Equation  (4.7),  we  see  that  we  are  integrating  Cl  from  the  ground  up.  We  expect  that 
the  turbulence  near  the  ground  is  the  major  contributor  to  the  phase  fluctuations.  On  the  other  hand. 
Equation  (4.6)  weights  the  integral  contributions  by  the  square  of  the  altitude.  Hence  we  expect  the 
higher  layers  of  turbulence  to  make  the  major  contribution  to  amplitude  fluctuation. 

The  integral-weighting  situation  should  tend  to  reduce  the  correlation  between  amplitude  and 
phase  fluctuations.  In  other  words,  the  instantaneous  mathematical  correlation  of  these  fluctuations 
should  be  low  because  they  arise  from  largely  different  physical  sources  of  turbulence,  different  air 
masses.  Incidentally,  Fried^  argues  that  the  difference  processes  for  phase  and  amplitude  are  uncorre¬ 
lated  due  to  isotropy.  Since  they  are  Gaussian,  it  follows  that  they  are  independent.  This  is  not  exactly 
the  same  point  we  are  making  here. 
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It  is  quite  possible  for  amplitude  and  phase  fluctuations  to  be  highly  correlated  in  the  rms  (average) 
sense,  e.  g.,  periods  of  large  phase  fluctuation  corresponding  to  large  amplitude  fluctuation.  The  issue  is 
whether  or  not  the  turbulent  activity  in  the  lower  layers  generally  tracks  the  turbulence  in  the  upper  lay¬ 
ers.  If  so,  devices  which  measure  amplitude  fluctuations  such  as  stellar  scintillation  should  be  useful  for 
characterizing  phase  fluctuations  most  of  the  time.  If  not,  we  may  notice  discrepancies  when  we 
characterize  the  atmosphere  with  a  scintillometer  in  an  effort  to  predict  the  performance  of  devices 
which  correct  phase. 

4.4.2  Fried’s  Parameter 

In  1964  Hufnagel  and  Stanley^  built  on  Tatarski’s  foundation  by  considering  the  long-term  average 
MTF  of  turbulence.  Their  method  involved  calculating  average  mutual  coherence  functions.  Fried^  fol¬ 
lowed  with  a  more  important  paper  calculating  both  long-term  and  short-term  average  MTFs.  He 
showed  that  there  is  much  less  high  spatial  frequency  attenuation  in  the  short-term  case.  In  addition,  he 
introduced  a  new  parameter  for  characterization  of  turbulence. 

Briefly,  the  mutual  coherence  function  of  the  wavefront  after  passing  through  turbulence  is  central 
to  the  analysis.  This  turns  out  to  be  what  Fried  calls  the  wave-structure  function 

=  +  Ds(f).  (4.8) 

This  quantity  had  previously  popped  up  in  Tatarski’s  work  in  the  derivation  of  Ds(i^  and  is  given  by 

D(f)  =  Al^'\  (4.9) 


Fried  defined  a  new  constant 

3/5 

(4.10) 


which  lets  us  write 

D(l^  =  6.88 

The  advantage  of  characterizing  turbulence  using  instead  of  A  is  the  physical  picture  in  which  turbu¬ 
lence  converts  the  effective  resolution  of  a  large  telescope  into  the  resolution  of  a  much  smaller  tele¬ 
scope.  is  the  diameter  of  that  smaller  telescope,  e.  g.,  a  diffraction-limited  telescope  of  diameter 
gives  resolution  comparable  to  a  very  large  telescope  through  the  given  turbulence. 

Since  Fried’s  work,  I'o  has  become  a  standard  measure  of  atmospheric  turbulence  or  quality  of 
astronomical  seeing.  It  is  a  useful  and  easily  understood  single  quantity,  but  of  course  it  falls  short  of 
completely  describing  the  three-dimensional  atmosphere.  One  should  also  keep  in  mind  that  the  use  of 
Kg  assumes  isoplanicity,  and  attempts  to  image  over  excessive  angles  will  violate  this  condition. 

4.5  SUMMARY 

When  Silverman  translated  Tatarski’s  book  in  1961 ,  a  major  contribution  to  the  mathematical  anal¬ 
ysis  of  turbulence  became  accessible  to  the  English  speaking  world.  This  approach  was  rapidly  adopted 


5/3 

(4.11) 


20 


by  many  investigators,  and  optical  propagation  through  turbulence  achieved  a  much  firmer  mathematical 
basis,  although  the  mathematical  difficulties  involved  in  exactly  calculating  some  of  the  functions  of 
interest  were,  and  are,  formidable. 

Fried  showed  that  the  high  spatial  frequency  transmission  through  turbulence  was  greatly 
attenuated  by  taking  long  term  averages  of  received  intensity,  which  explained  some  anomalous  results 
of  imaging  experiments.  In  the  short  term  case,  there  were  different  effects  for  the  near-field  and  far- 
field  cases,  but  in  the  long  term  case  this  distinction  did  not  exist.  Accordingly,  Fried  suggested  that 
propagation  be  studied  using  long  exposures  so  that  the  results  could  be  better  understood.  He  did  hold 
out  the  hope  that  short  exposures  could  be  more  useful  after  theory  advanced. 

Labeyrie^  was  the  first  to  show  how  short  exposures  could  be  used  to  extract  information  from  the 
higher  spatial  frequencies.  The  experimental  protocols  which  followed  were  able  to  sidestep  theoretical 
understanding  of  turbulence  somewhat.  Long-term  averages  of  quantities  derived  from  short  exposures 
tended  to  reduce  the  need  for  understanding  instantaneous  turbulence.  Measuring  the  effect  of  turbu¬ 
lence  by  observing  stars  enabled  corrections  to  be  made  without  any  understanding  of  the  underlying 
phenomenology.  Experimentalists  plunged  into  areas  which  were  thought  on  theoretical  grounds  to  be 
quicksand,  sometimes  surprising  the  experts. 

We  will  consider  these  issues  in  detail  below,  but  first  we  will  consider  the  interesting  and  rela¬ 
tively  hardware-intensive  subject  of  adaptive  optics.  While  speckle  imaging  is  perhaps  superior  to  adap¬ 
tive  optics  for  dim  astronomical  objects,  we  will  find  that  some  earth-orbiting  satellites  may  be  more 
easily  imaged  adaptively.  Considering  speckle  imaging  and  adaptive  optics  together  will  suggest  inter¬ 
esting  combinations  of  these  techniques. 
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5.  ADAPTIVE  OPTICS 


5.1  OVERVIEW 

Adaptive  or  active  optics  describes  any  system  in  which  the  optical  components  change 
configuration  during  operation.  Historically,  due  to  the  small  wavelengths  involved,  most  man-made 
optical  devices  were  rigid  lenses  and  mirrors  made  out  of  a  material  like  glass.  Such  devices  operate  in 
an  “open  loop”  manner  with  fixed  configuration  unaffected  by  equipment  operation. 

Biological  imaging  systems  have  not  been  bound  by  such  constraints;  the  human  eye  is  able  to 
adjust  its  focus  by  changing  the  shape  of  the  lens,  maintaining  focus  on  moving  objects.  Another  famil¬ 
iar  example  of  adaptive  optics  is  a  camera  with  automatic  focusing.  Here  the  adaptation  is  performed  by 
moving  rigid  lenses.  In  both  cases  a  feedback  control  system  operates  to  minimize  deviation  from 
desired  performance. 

There  are  other  examples  of  applications  in  which  relatively  slow  adaptation  is  effective,  and  in 
which  the  degradation  can  be  characterized  by  a  small  number  of  parameters.  Large  space-based  optics 
will  require  some  adaptation  to  maintain  their  configuration  in  the  presence  of  nonuniform  solar  heating. 
Large  land-based  telescopes^^  constructed  of  many  smaller  mirrors  will  use  adaptive  methods  to  reduce 
required  rigidity  and  weight  while  decreasing  optical  errors. 

The  most  challenging  case  is  adapting  to  rapidly  varying  optical  disturbances  which  are  them¬ 
selves  complex  in  nature,  as  in  the  important  problem  of  atmospheric  turbulence.  This  was  first 
achieved  about  15  years  ago  and  is  our  primary  concern.  We  will  henceforth  use  the  terminology  adap¬ 
tive  optics  to  refer  to  such  systems.  Hardy  ^  has  written  an  excellent  summary  of  the  area  as  of  about 
1978.  Most  of  it  is  still  highly  relevant,  although  of  course  there  have  been  certain  advances  since  that 
time. 

We  can  divide  adaptive  optics  systems  into  two  categories,  active  and  passive.  The  active  is  con¬ 
cerned  with  transmitting  energy,  typically  produced  by  a  laser,  through  turbulence.  The  beam  is  pre¬ 
compensated  to  largely  cancel  the  effects  of  the  turbulence.  The  passive  system  corrects  the  received 
wavefront  to  cancel  turbulence.  We  will  say  no  more  about  active  systems. 

There  are  two  basic  approaches  to  passive  adaptive  optics,  illustrated  in  Figure  5.1.  In  the 
front  compensation  approach,  a  wavefront  modifier  attempts  to  compensate  for  the  effects  of  atmos¬ 
pheric  turbulence.  After  modification  by  this  device,  the  wavefront  is  split  so  that  part  of  the  energy  is 
used  to  make  an  image  and  the  rest  is  sent  to  a  wavefront  detector.  The  latter  device  maps  out  the  resid¬ 
ual  wavefront  deviation,  which  is  then  converted  to  modifier  control  signals  by  the  processor.  These  sig¬ 
nals  are  added  to  the  existing  ones  and  the  closed-loop  system  attempts  to  null  the  error. 

In  the  image  sharpening  approach  deliberate  distortions  are  introduced  in  the  received  wavefront 
over  a  wide  variety  of  amplitudes  and  shapes.  The  intensity  detector  determines  which  of  the  trial 
modifications  introduced  by  the  aperture  tagger  is  effective  in  improving  image  “sharpness.”  The  pro¬ 
cessor  converts  that  information  into  control  signals  which  add  an  appropriate  modification  to  the  wave- 
front  through  the  wavefront  modifier. 

The  most  successful  approach  to  date  has  been  the  wavefront  compensation  one,  which  we  will 
henceforth  refer  to  as  compensated  imagings  and  we  will  spend  most  of  our  time  discussing  it.  The 
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(a)  Wavefront  Compensation 


(b)  Image  Sharpening 


Figure  5.1.  Adaptive  optics  systems. 


image  sharpening  work  is  important  in  some  applications  such  as  solar  imaging  where  compensated 
imaging  is  inappropriate.  More  importantly,  the  reported  work  on  image  sharpening  will  influence  our 
analysis  of  compensated  image  restoration.  It  will  also  be  of  interest  in  our  discussion  of  speckle  imag¬ 
ing  of  satellites. 

5.2  HISTORY 

In  the  early  1970s  DARPA  funded  experimental  investigations  into  compensated  imaging  by 
Hardy  and  his  colleagues  at  Itek.  In  late  1973  a  practical  compensated  imaging  device^  was  demon¬ 
strated.  We  will  briefly  summarize  the  configuration  of  this  device  here;  we  will  describe  the  compo¬ 
nents  in  more  detail  below. 

A  30  cm  aperture  was  divided  into  21  subapertures.  The  wavefront  modifying  device  was  a  mono¬ 
lithic  piezoelectric  mirror  (MPM).  The  MPM  itself  was  not  as  large  as  30  cm;  the  light  was  collimated 
by  a  lens  before  striking  the  smaller  MPM.  The  wavefront  detector  was  an  array  of  interferometers, 
which  produced  32  wavefront  tilt  measurements,  16  in  each  of  two  orthogonal  directions.  An  analog 
computer  solved  the  least  squares  problem  of  converting  the  32  available  measurements  into  the  best 
deflection  for  each  of  the  21  MPM  actuators.  In  addition,  a  separate  system  processed  the  wavefront 
sensing  data  and  corrected  overall  waveform  tilt.  This  was  a  great  step  forward,  and  required  advances 
in  almost  every  component  to  work  successfully. 
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The  next  milestone  was  a  69  subaperture  MPM-based  system  built  by  Hardy  at  Itek  Corporation 
for  Division  5  of  Lincoln  Laboratory.  A  60  cm  aperture  was  used,  so  the  subapertures  were  slightly 
larger  than  in  Hardy’s  original  design.  The  same  type  of  interferometer  was  used  for  phase  sensing, 
although  its  wavelength  insensitivity  was  not  needed  for  this  laser  application.  This  appears  to  be  very 
much  a  scaled*up  version  of  its  predecessor.  The  only  apparent  improvement  besides  the  increased  size 
is  the  replacement  of  the  analog  phase-lock-loop  phase  measurement  circuitry  by  digital  correlation 
equipment. 

The  next  step  in  the  progression  is  an  imaging  system known  as  the  Compensated  Imaging  Sys¬ 
tem  (CIS),  which  was  put  into  service  in  Hawaii  in  1981.  This  device  was  again  built  by  Itek,  using 
essentially  the  same  technology  as  its  predecessors,  but  with  a  1.6  m  aperture.  The  CIS  is  operated  by 
AVCO  as  part  of  AMOS/MOTIF  (ARPA  Maui  Optical  Station/Maui  Optical  Tracking  and  Identification 
Facility).  This  is  a  large  and  impressive  piece  of  equipment;  it  is  currently  the  logical  choice  for  experi¬ 
ments  in  adaptive  optical  space  surveillance. 

The  astronomy  community  has  not  been  able  to  afford  compensated  imaging  equipment  for  the 
most  part.  The  expense  and  the  limitation  to  bright  objects  has  made  the  payoff  look  less  than  the  pay¬ 
out.  There  is  a  potential  program*^  of  interest  being  promoted  by  the  NO  AO  organization.  This  would 
perform  adaptive  correction  in  the  infrared,  the  reasoning  being  that  a  large  aperture  requires  many  fewer 
subapertures  due  to  the  increase  in  in  the  infrared. 

The  data  which  might  be  produced  by  such  a  system  are  not  of  as  much  interest  as  the  theoretical 
papers  written  by  the  university  astronomers  involved  in  the  project.  It  seems  that  when  you  are  knee- 
deep  in  data  from  a  real  system,  you  spend  most  of  your  time  looking  at  lists  of  satellites  and  sending 
videotapes  to  sponsors.  The  astronomers,  on  the  other  hand,  lacking  data  from  any  operational  system, 
are  forced  to  analyze  their  projected  system  in  meticulous  detail  to  minimize  risk  and  cost.  Where  the¬ 
ory  fails,  they  usually  have  access  to  computers  suitable  for  simulation.  We  will  reference  some  of  their 
work  when  we  discuss  errors  in  adaptive  systems. 

Departing  from  the  mainstream  of  compensated  imaging,  Muller  and  Buffington  published  some 
results*^  on  image  sharpening  at  the  same  time  that  Itek  was  working  on  the  initial  prototype  of  their 
equipment.  They  showed  that  the  integral  of  the  squared  image  intensity  is  a  suitable  sharpness  measure 
for  controlling  adaptive  optics.  They  considered  other  measures  as  well,  reporting  the  results  of  com¬ 
puter  simulations.  Later  Buffington  et  described  experiments  with  real  hardware. 

5.3  MPM  DESCRIPTION 

The  monolithic  piezoelectric  mirror  (MPM)  is  the  most  successful  wavefront  modifying  device  to 
date.  Hardy  et  al^  describe  a  21 -element  MPM  in  detail.  This  device  consists  of  a  monolithic  cylinder 
of  sintered  lead-zirconate-titanate  (PZT),  a  standard  piezoelectric  material.  The  dimensions  are  about 
5  cm  diameter  and  1.2  cm  thick. 

An  array  of  electrodes  is  placed  on  the  top  of  the  PZT,  addressed  by  wires  passing  through  to  the 
back.  A  common  electrode  covers  the  back  surface  and  a  mirror  is  glued  to  the  front  surface. 

When  a  positive  or  negative  voltage  is  applied  to  an  electrode,  a  small  localized  region  of  the  mir¬ 
ror  deflects  in  the  corresponding  direction.  There  is  reported  to  be  little  interaction  between  electrodes,  a 
major  advantage  of  this  type  of  wavefront  corrector. 
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One  disadvantage  is  the  high  voltage  required  to  operate  the  MPM;  2-4  kV  to  deflect  the  required 
2000  nm  or  so.  The  inherent  frequency  response  of  the  material  is  from  dc  to  over  10  kHz,  but  the 
arrangement  of  amplifiers  used  in  this  original  compensated  imaging  system  reduced  that  to  about 
1  kHz.  This  does  not  represent  any  technology  limitation  today. 

5.4  WAVEFRONT  SENSOR  DESCRIPTION 

One  of  the  problems  with  compensated  imaging  is  that  it  requires  enough  light  in  each  subaperture 
to  estimate  wavefront  distortion.  If  the  interferometers  required  filters  to  produce  narrowband  light,  the 
situation  would  be  exacerbated.  Hence  practical  interferometers  for  this  application  operate  with  white 
light. 

If  we  imagine  the  apparatus  required  for  Young’s  experiment,  we  will  be  totally  unprepared  for  the 
variety  and  complexity  of  modem  interferometers.  Rather  than  delve  into  the  general  subject,  we  will 
confine  our  discussion  to  the  ac  shearing  interferometers  used  in  the  series  of  Itek  interferometers  and 
described  by  Hardy^’  ’  and  Koliopoulos.'^ 

First  we  will  be  sure  we  understand  just  what  we  mean  by  wavefront  “tilt.”  We  assume  that  the 
effects  of  turbulence  are  continuous,  so  that  optical  path  lengths  do  not  vary  significantly  for  very  small 
spatial  separations.  In  particular,  the  spatial  path  length  function  could  be  expanded  in  some  small 
subaperture  in  a  power  series  and  represented  by  some  constant  plus  the  linear  term,  or  actually  one  lin¬ 
ear  term  for  the  x  direction  and  one  for  y. 

Now  if  we  separate  two  pinholes  in  a  Young’s  experiment  by  some  small  A  (a  subaperture  size), 
we  will  get  an  interference  pattern  which  tells  us  the  path  difference,  or  equivalently  the  waveform  tilt. 
A  two-dimensional  array  of  x  and  y  tilts  is  really  an  expression  of  the  differential  path  length  function  in 
terms  of  multiple  truncated  power  series.  If  the  subapertures  are  small  enough,  this  is  a  highly  accurate 
representation. 

In  fact.  Young’s  experiment  would  tell  us  the  phase  difference  of  the  incident  narrowband  light 
rather  than  the  length  difference.  However,  since  the  optical  effects  of  turbulence  are  caused  by  varia¬ 
tions  in  refractive  index,  the  cause  of  the  phase  difference  may  be  viewed  as  an  underlying  constant  path 
length  difference.  Hence  the  phase  difference  will  be  proportional  to  frequency  and  white  light  inter¬ 
ferometers  can  measure  the  distance,  or,  as  in  our  case,  the  wavefront  slope  over  the  subaperture. 

5.4.1  Hardy’s  Interferometer 

We  will  now  describe  Hardy’s  variable-shear  rotating  ac  interferometer,  shown  in  Figure  5.2, 
which  is  about  as  complicated  as  its  name.  This  device  was  developed  especially  for  compensated  imag¬ 
ing;  in  fact  the  first  one  was  used  in  an  upgrade  of  the  original  compensated  imaging  equipment. 

Typically  in  a  compensated  imaging  system,  the  light  from  the  full  aperture  is  collected  and  col¬ 
limated  into  a  much  smaller  diameter  suitable  for  modification  by  the  MPM.  After  MPM  correction  but 
before  entering  the  final  lens  to  be  focused  on  the  imaging  system,  the  light  is  split  into  two  paths,  one 
for  the  X  tilt  measurement  and  one  for  the  y  tilt.  Lens  LI  focuses  the  light  onto  G,  a  rotating  circular 
grating  with  alternating  clear  and  opaque  radial  lines  of  equal  thickness.  A  point  of  light  projected  any¬ 
where  on  this  grating  produces  a  square  wave  time  function  on  the  other  side.  After  diffraction  by  the 
grating,  lens  L2  forms  multiple  images  of  the  wavefront  on  the  detector  array.  When  the  aperture  size 
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(and  number  of  subapertures)  is  increased,  the  detector  array  and  associated  electronics  are  the  only 
things  that  need  to  be  modified. 

We  can  view  the  lens  arrangement  as  an  optical  Fourier  transformer;  when  we  multiply  the  Fourier 
transform  of  the  wavefront  by  a  square  wave  that  is  the  same  as  convolving  the  wavefront  by  an  impulse 
train  of  odd  harmonics.  Not  shown  in  Figure  5.2  is  what  happens  after  detection;  the  detector  signals  are 
processed  by  electrical  filters  which  remove  unwanted  harmonics.  Even  though  only  three  images  are 
shown  in  Figure  5.2,  there  are  practically  an  infinite  number  of  them.  Of  course  for  a  given  shear  (angu¬ 
lar  separation)  only  a  finite  number  of  them  actually  overlap  in  the  area  of  the  detector.  The  electrical 
filtering  is  necessary  to  sort  them  out  so  that  a  clean  measurement  may  be  made. 

If  you  view  the  gratings  from  the  front,  the  reference  source  is  focused  at  say  12  o’clock  on  Gx  and 
3  o’clock  on  Gy.  This  gives  the  shear  in  the  right  direction.  The  shear  is  varied  by  moving  Gx  left  and 
right  and  Gy  up  and  down  so  that  the  spatial  frequency  of  the  square  wave  diffraction  pattern  decreases 
or  increases.  This  has  the  effect  of  varying  the  gain  in  the  feedback  control  system,  and  helps  compen¬ 
sate  for  various  target  and  atmosphere  conditions. 

5-4.2  Wavelength  Dependence 

Now  that  we  understand  schematically  how  it  works,  let  us  consider  the  white  light  property.  As 
Hardy*  explains,  the  shear  0^  between  the  overlapping  images  is  defined  by  the  equation 

sin05  =±nXv  /2=1,3,5,...  (5.1) 

where  \  is  the  optical  wavelength  and  v  is  the  grating  spatial  frequency. 

Denote  the  slope  of  the  incoming  wavefront  a(jc),  and  the  focal  length  of  lens  L2  F.  The  shear  dis¬ 
tance  in  the  detector  plane  is  given  by  5  =  0^^  and  we  have  the  following  expression  for  the  phase 
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(5.2) 


difference  between  the  sheared  wavefronts  in  the  detector  plane: 

(t)(A:)  =  a(A)|-. 

The  interference  phase  angle  is  the  quantity  of  interest  when  the  multiple  images  interfere.  This 
angle  is  invariant  to  X  since  S  is  proportional  to  X  for  small  angles.  Hence  is  proportional  to  the 
wavefront  slope,  a  dimensionless  quantity,  and  the  system  operates  in  white  light. 

5.4.3  Electronics 

Each  detector  outputs  an  ac  signal  in  the  kHz  region  whose  phase  is  proportional  to  wavefront  tilt 
in  the  corresponding  subaperture.  We  have  already  mentioned  that  the  extraneous  sidebands  are  filtered 
out  at  this  point.  The  remaining  processing  is  amplification  and  AGC  followed  by  phase  detection  with 
a  phase-locked-loop.  X  and  Y  phase  references  are  generated  from  the  rotating  gratings  for  use  in  this 
process. 

For  further  details  on  this  device,  consult  the  references.  Note  that  while  some  interferometers 
measure  path  lengths  directly,  this  one  measures  wavefront  slope  instead.  That  is  no  major  problem,  but 
it  certainly  influences  the  processor,  discussed  next. 

5.5  PROCESSOR 

The  detector  array  in  a  compensated  imaging  system  typically  is  a  square  grid  inside  a  circular 
aperture,  as  shown  in  Figure  5.3.  The  dots  represent  the  MPM  actuators;  the  actual  physical  locations  in 
the  image  of  the  x  and  y  tilt  measurements  are  between  the  MPM  actuators  about  where  the  As  are  in  the 
figure.  We  have  x  and  y  phase  differences,  so  what  we  would  like  to  do  is  perform  a  least-squares  fit  of 
the  data  to  solve  this  overdetermined  system  of  equations  for  the  best  actuator  displacements. 

Hudgin^^  addressed  this  issue,  concluding  that  both  analog  and  digital  processors  were  feasible. 
When  one  reads  his  paper,  one  wonders  why  he  traveled  to  Fourier  space  rather  than  just  solving  the  nor¬ 
mal  equations  directly.  Herrmann^^  evidently  had  the  same  question,  and  he  provided  a  more  coherent 
analysis  of  the  problem.  He  also  corrected  a  couple  of  Hudgin’s  errors.  Hunt^^  also  gave  a  general  treat¬ 
ment  of  the  problem  in  an  explosion  of  matrix  equations. 

Hudgin  used  certain  approximations  valid  for  large  arrays  in  his  derivation,  so  it  is  not  guaranteed 
optimum  for  finite  arrays.  His  solution  for  an  unknown  phase  involves  the  average  of  its  four  nearest 
neighbors  plus  a  weighted  sum  of  the  four  connecting  phase  differences. 

The  analog  realization  consists  of  a  resistive  network  and  many  operational  amplifiers;  this  is  the 
approach  taken  in  the  current  compensated  imaging  systems.  The  digital  approach  involves  adders  only; 
divisions  are  performed  by  shifting. 

The  exact  solution  is  a  least  squares  matrix  problem.  We  can  define  a  constant  connection  matrix 
A,  relating  the  phase  difference  vector  A  to  the  phase  vector  (])  by  A  A  =  (j).  This  matrix  A  consists  of 
mostly  zeros  with  a  1  and  a  -1  in  each  row. 

The  solution  to  this  least  squares  problem  is  well  known;  the  inverse  of  the  “sample  covariance 
matrix”  A 'A  is  multiplied  by  A 'A  in  real  time  to  obtain  the  appropriate  phases  (MPM  actuator  signals). 
Since  A  doesn’t  change  [A  'A]"'  need  only  be  precomputed  once. 
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Figure  5  J.  Wavefront  reconstruction. 

The  analog  version  requires  additional  feedback  loops  to  ensure  that  the  tilt  is  entirely  applied  to 
the  tilt  equipment  rather  than  using  up  the  travel  distance  of  the  MPM.  Similarly,  the  average  deflection 
of  the  MPM  is  kept  near  zero,  since  absolute  path  distances  are  not  important.  These  factors  are  more 
easily  accommodated  using  digital  techniques. 

There  are  many  issues  which  can  be  brought  up  such  as  noise  propagation  and  time  delay,  but  they 
would  take  us  too  far  afield.  Current  designs  mostly  use  analog  processors,  but  it  seems  likely  that 
future  systems  will  employ  digital  processors,  probably  employing  the  optimum  algorithm  rather  than 
Hudgin’s  approximations. 
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6.  SPECKLE  IMAGING 


6.1  SPECKLE  INTERFEROMETRY 

6.1.1  Labcyrie’s  Method 

In  1970  Labeyrie^  published  a  short  paper  in  which  he  pointed  out  that  the  diffraction-limited 
power  spectrum  of  the  intensity  of  an  astronomical  object  could  be  estimated  using  a  sequence  of  short- 
exposure  images.  Fried^  had  already  showed  that  short-exposure  images  had  modulation  transfer  func¬ 
tions  (MTFs)  which  were  significant  out  to  the  diffraction  limit.  He  did  not,  however,  give  any  algo¬ 
rithm  for  converting  a  short  term  image  or  image  sequence  into  a  practical  high-resolution  image. 

Labeyrie  did  not  solve  the  high-resolution  imaging  problem  either,  but  his  insight  was  that  rather 
than  dealing  with  the  standard  incoherent  quasi-monochromatic  imaging  equation  Equation  (3.32)  one 
could  instead  profitably  work  with  its  magnitude  squared 

|/(«,v)|-  =  |0(m,v)|2  \t{u,v)\^.  (6.1) 

Recall  that  T  (zv,v)  is  the  transfer  function  of  the  atmosphere  plus  lens  combination.  Since  there  is  a  ran¬ 
dom  phase  added  by  turbulence,  the  mean  value  of  averages  of  /(w,v)  tends  to  zero  for  spatial  frequen¬ 
cies  beyond  1  cycle/arc  sec  or  so.  If  one  considers  instead  the  magnitude  squared  of  I{u,v),  the  mean 
value  remains  significant. 


6.1. 1.1  Frequency  Domain  Intuitive  Explanation 


Before  discussing  this  more  rigorously,  consider  the  heuristic  analysis  of  Roddier.^^  Recall  our  dis¬ 
cussion  of  the  interferometric  view  of  imaging.  An  image  can  be  imagined  to  be  made  up  of  the  super¬ 
position  of  fringe  patterns  caused  by  pairwise  combinations  of  points  in  the  exit  pupil.  In  the  same  way, 
we  can  imagine  an  aperture  to  be  made  up  of  subapertures  larger  than  a  point;  the  whole  image  can  be 
considered  to  be  a  superposition  of  fringes  from  the  pairwise  combinations  of  subapertures. 

As  a  point  of  departure,  again  consider  the  form  of  T (w,v)  through  a  diffraction-limited  lens.  As 
we  found  in  Section  3.6,  defining  f  to  be  our  vector  frequency  (w,  v). 


nD= 


N{0) 


(6.2) 


where  N  (f)  the  pupil  redundancy  is  the  number  of  subaperture  pairs  which  add  to  produce  frequency  f. 
This  is  again  proportional  to  the  shaded  area  in  Figure  3.5. 

Now  consider  the  situation  of  imaging  through  the  atmosphere.  Let  our  subaperture  size  be  and 
let  the  random  phase  delay  6^  on  subaperture  k  be  constant  and  independent  of  the  phase  delay  on  the 
other  subapertures.  We  can  write  our  imaging  equation  as 


/(f)  =  0(0  5(0 

where  the  transfer  function  S(0  is  given  by  the  random  sum 

5(0  = 


1 


y  € 
A'(O) 


(6.3) 


(6.4) 
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The  sum  is  over  all  subaperture  pairs  contributing  to  frequency  f  and  the  phase  angle  is  the  difference  in 
phase  between  each  such  pair  of  subapertures. 

The  situation  is  mathematically  a  random  walk.  The  average  magnitude  I  S(f)  I  grows  as  the 
square  root  of  the  number  of  pairs.  We  conclude  that  the  average  power  spectrum  (  I  5(0  I  is  propor¬ 
tional  to  A^(0  rather  than  A^^(0. 

This  analysis  gives  us  helpful  insight,  but  unfortunately  in  practice  speckle  imaging  extracts  these 
frequencies  through  the  intermediate  step  of  periodograms.  The  poor  statistical  properties  of  these  esti¬ 
mators  produce  a  large  variance  which  worsens  things.  We  will  discuss  this  situation  in  more  detail 
below. 


6.1. 1.2  Spatial  Domain  Intuitive  Explanation 

There  is  also  a  very  useful  spatial  domain  heuristic  analysis  which  we  will  discuss  in  this  section. 
It  originates  in  the  classical  astronomy  problem^^  of  identifying  and  analyzing  binary  st£u*s.  Assume  we 
have  an  image  of  two  objects  separated  by  distance  d  such  as  depicted  in  Figure  6.1a.  The  point  spread 
function  of  typical  turbulence  tends  to  produce  multiple  copies  of  the  object  being  viewed  at  random 
positions  within  the  fuzzy  disk  (Airy  disk)  corresponding  to  the  long-term  exposure.  This  is  shown  in 
Figure  6.1b.  This  effect  is  quite  evident  if  bright  stars  are  photographed  properly  with  a  large  instrument 
and  a  short  exposure. 


(a)  Diffraction-limited  Image 


(b)  Short  Exposure  Speckle  Image 


Figure  6.1.  Schematic  speckle  imaging. 


Now  rather  than  consider  the  effect  of  the  turbulence  in  the  frequency  domain,  in  terms  of  MTFs 
and  periodograms,  consider  the  equivalent  autocorrelation  function.  The  (ideal)  autocorrelation  function 
of  the  image  in  Figure  6.1a  consists  of  the  superposition  of  the  autocorrelation  functions  of  the  two  iso¬ 
lated  objects  at  the  origin  and  two  cross-correlation  terms  at  shifts  of  d. 

Now  consider  the  autocorrelation  function  of  Figure  6.1b.  It  is  obvious  that  for  small  shifts  on  the 
order  of  the  diameter  of  the  individual  objects  we  will  get  a  pretty  good  estimate  of  the  ideal  autocorrela¬ 
tion.  The  only  errors  will  occur  when  the  multiple  images  are  spaced  less  than  one  of  their  diameters 
apart.  Similarly,  the  autocorrelation  for  shifts  of  d  will  be  quite  accurate. 
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The  main  errors  occur  for  other  shifts  when  spurious  cross-correlation  terms  appear.  Recalling  that 
the  positions  of  the  images  of  Figure  6.1b  are  random,  it  is  clear  that  averaging  many  of  these  autocorre¬ 
lation  functions  will  tend  to  reduce  the  amplitude  of  these  spurious  correlations.  Hence  we  obtain  the 
proper  diffraction-limited  autocorrelation  function  as  Labeyrie  proposed.  This  argument  is  of  course 
more  convincing  for  isolated  small  objects  than  extended  objects. 

6.L1.3  Summary 

Now  let  us  summarize  the  steps  needed  in  Labeyrie’s  method.  Note  that  the  original  algorithm 
used  optical  processing  due  to  the  unavailability  of  suitable  digital  computers.  We  will  accordingly 
update  the  method  for  modem  equipment. 

We  first  obtain  a  short  exposure  image  of  the  object  under  study.  We  compute  the  two-dimensional 
discrete  Fourier  transform  (DFT).  We  then  take  the  square  of  this  DFT  as  a  single  estimate  of  the  image 
power  spectrum.  This  estimator,  called  a  periodogram,  is  known^"^  to  be  a  very  poor  one  in  terms  of  bias 
and  variance.  While  it  is  possible  to  improve  the  estimator  by  certain  smoothing  processes,  that  is  obvi¬ 
ously  not  what  we  desire  since  high  resolution  is  our  goal.  The  practical  solution  is  to  average  many 
periodograms,  as  Labeyrie  proposed. 

There  is  one  final  step  in  the  method.  Since  the  atmosphere/lens  combination  has  modified  the 
spatial  frequencies  by  some  MTF,  to  restore  the  object  power  spectrum  (autocorrelation)  it  is  necessary 
to  make  a  correction,  as  discussed  below.  This  will  be  a  critical  step  in  later  imaging  algorithms  and  in 
some  cases  the  most  difficult  step  to  perform. 

6.1.2  OTF  Corrections 

6.1.2.1  Theoretical  Analysis 

Now  that  we  have  a  heuristic  picture  of  the  turbulence  MTF,  we  will  consider  briefly  the  rigorous 
analysis.  Recall  that  in  1966  Fried^  analyzed  both  the  long  term  and  short  term  imaging  MTFs.  Mathe¬ 
matically,  he  evaluated  (S(f))-  Four  years  later  Labeyrie  published  his  paper,  and  there  was  a  need  to 
evaluate  (  I  5(f)  1  ^);  this  was  done  by  Korff*^  in  1973. 

This  was  a  rather  difficult  calculation  involving  the  fourth  order  statistics  of  the  amplitude  of  5. 
While  certain  approximations  were  made,  the  results  have  been  shown  to  agree  quite  well  with  experi¬ 
ment.  We  will  not  delve  into  the  mathematics  here,  although  it  is  not  too  hard  to  follow  if  you  are  famil¬ 
iar  with  Tatarski’s  book  and  Fried’s  paper.  The  bottom  line  is  not  even  an  equation,  since  Korff  was 
unable  to  solve  the  problem  analytically,  but  a  graph  reflecting  the  results  of  computer  calculations. 

Figure  6.2  reproduces  Korff  s  results.  We  can  see  the  substantial  improvement  in  high  spatial  fre¬ 
quency  response  using  Labeyrie’s  method.  Notice  particularly  that  for  the  large  diameter  case,  which 
might  correspond  to  a  4  m  telescope,  the  Labeyrie  transfer  function  seems  likely  to  support  at  least  an 
order  of  magnitude  better  resolution.  The  exact  figure  depends  on  signal-to-noise  considerations  and  so 
will  differ  from  object  to  object.  Of  course  as  the  seeing  varies,  even  with  constant  telescope  diameter, 
different  transfer  function  curves  are  appropriate. 

In  the  case  of  small  diameters,  it  appears  from  Figure  6.2  that  a  MTF  correction  might  improve 
long-term  exposure  data  as  well.  Goodman  and  Belsher^^  analyzed  the  situation  and  concluded  that  this 
was  possible. 
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Figure  6.2.  Turbulence  MTF  for  Labeyrie  (solid),  long-term  exposure  (dotted). 

6.1.2.2  Using  Reference  Stars 

If  there  is  one  word  which  describes  the  theory  of  speckle  imaging  it  is  probably  “approxima¬ 
tion.”  From  the  first  principles  of  optics  through  the  Van  Cittert-Zemike  theorem  and  the  models  of  tur¬ 
bulence  we  have  seen  terms  such  as  quasi-monochromatic,  incoherent,  plane  wave,  isoplanatic,  Kolmo¬ 
gorov  model,  etc.  Many  of  these  approximations  are  quite  good  in  the  proper  circumstances  but  there  is 
a  gradual  erosion  of  confidence  in  the  theory  as  they  pile  up  one  after  the  other.  We  still  feel  that  our 
intuitive  understanding  is  quite  well  served  by  these  concepts,  but  it  is  logical  to  wonder  how  far  a  real 
system  transfer  function  deviates  from  the  theory. 

There  are  also  cases  in  which  some  of  the  required  assumptions  are  not  met,  and  we  have  no  reason 
to  believe  that  the  theory  is  applicable.  We  might  still  have  data  available  and  want  to  restore  its  resolu¬ 
tion  as  far  as  possible.  What  is  needed  is  some  metric  which  can  guide  us  in  our  correction  for  the 
unknown  MTF  of  the  total  imaging  process. 

This  situation  has  been  recognized  from  the  beginning,  and  was  even  more  important  then  since 
Korff  s  paper  and  others  had  not  yet  been  published.  One  solution,  mentioned  briefly  by  Labeyrie,^  is  to 
experimentally  measure  the  long-term  average  of  the  short-term  MTF  of  the  atmosphere  plus  imaging 
system.  This  is  conveniently  done  using  an  unresolvable  star,  since  such  a  point  object  produces  an 
image  which  is  the  point  spread  function  (PSF)  itself,  i.  e.,  the  Fourier  transform  of  the  OTF.  The  phases 
of  the  OTF  are  not  recoverable  due  to  the  turbulence,  but  the  average  magnitude  gives  us  a  useful  MTF 
correction. 

One  problem  occurs  when  the  nearest  suitable  reference  star  is  further  than  some  isoplanatic  angle 
from  the  object  under  study.  Theory  tells  us  that  the  PSFs  will  differ  in  this  case.  The  solution  lies  in 
long-term  averages.  While  a  short  sequence  of  degraded  objects  might  experience  a  different  MTF  than 
the  nearest  reference  star,  we  expect  that  a  long-term  average  will  differ  much  less.  Experiment  supports 
this  argument. 
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There  are  still  problems  with  reference  stars.  It  is  helpful  if  they  are  fairly  bright,  which  limits 
their  availability.  If  extremely  bright  objects  such  as  low-orbit  satellites  are  imaged,  they  may  require  a 
filter  to  avoid  saturation  which  requires  even  brighter  reference  stars.  One  very  interesting  approach^^  is 
to  create  “artificial  stars”  using  lasers.  Although  this  is  still  in  the  preliminary  investigation  stage,  there 
is  an  obvious  attraction  to  producing  “reference  stars  on  demand.” 

6.2  PHASE  RETRIEVAL  FOR  SPECKLE  IMAGING 

Up  to  now  we  have  discussed  speckle  interferometry,  in  which  we  cleverly  process  short -exposure 
images  to  obtain  near-diffraction-limited  autocorrelation  functions.  We  do  not,  however,  obtain  actual 
images.  As  Labeyrie  put  it,^  “Losing  the  phase  makes  it  impossible  to  reconstruct  the  object,  except  if  it 
has  a  center  of  symmetry.” 

Fortunately,  he  was  wrong  in  this  assertion;  there  are  methods  by  which  the  phase  can  be  retrieved 
from  the  data.  There  are  also  purely  mathematical  methods  in  which  the  phases  are  retrieved  from  the 
magnitudes  by  the  use  of  iterative  procedures  assuming  certain  constraints.  We  will  discuss  those  algo¬ 
rithms  next. 

6.2.1  Knox-Thompson  Method 

The  first  practical  phase  retrieval  method  is  due  to  Knox  and  Thompson.^^  This  important  paper 
showed  that  indeed  phase  information  did  exist  in  speckle  data,  and  provided  an  algorithm  to  retrieve  it. 
A  one-dimensional  simulation  is  the  only  result  shown,  although  Knox^^  later  provided  two- 
dimensional  simulation  results. 

Practically  any  algorithm  can  be  analyzed  in  both  the  time  (space)  domain  and  the  frequency 
domain.  One  is  often  clearer  than  the  other,  although  which  one  sometimes  depends  on  the  observer. 
We  have  tried  to  present  both  approaches  where  possible  in  this  report.  Knox  and  Thompson  explain 
their  algorithm  in  the  space  domain  by  recourse  to  autocorrelation  functions.  It  is  clearer  and  fits  in 
much  better  with  our  development  to  discuss  it  in  the  frequency  domain. 

Recall  the  imaging  equation  Equation  (6.3),  and  the  frequency  domain  intuitive  explanation  of 
incoherent  imaging  through  the  atmosphere.  We  can  view  our  transfer  function  S  (f)  as  being  approxi¬ 
mately  made  up  of  subapertures  of  dimension  about  which  have  nearly  constant  phase  shift.  The  shift 
is  of  course  random  over  large  distances. 

When  we  compute  the  DFT  of  an  image,  i(x,y)  — »/(w,v),  we  obtain  a  two-dimensional  array  of 
spatial  frequencies.  While  the  absolute  phases  of  these  frequencies  ^(u,v)  are  the  phases  of  the  true 
image  plus  some  random  phases,  if  two  phases  (|)(w  i,vi)  =  (J)(f| )  and  (t)(w2»^2)  =  ^>(^2)  closely  spaced 
in  frequency,  the  random  phase  components  are  nearly  identical.  Hence  we  can  take  the  difference 
between  such  closely  spaced  phases  and  effectively  obtain  the  difference  between  the  true  image  phases. 
Closely  spaced  here  means  that  the  length  of  the  frequency  vector  I  G  -f2  I  =  Af  (in  cycles/rad)  is  less 
than  Kfj/X. 

As  the  mathematicians  would  say,  the  problem  is  solved  because  we  have  reduced  it  to  an  equiva¬ 
lent  problem  already  solved.  That  is  of  course  the  processor  used  in  compensated  imaging,  and  the  cal¬ 
culation  depicted  in  Figure  5.3.  The  situation  is  the  same,  with  measured  differences  being  combined  to 
obtain  the  actual  function.  One  distinction  is  that  we  have  exactly  the  measurements  shown  in 
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Figure  5.3  in  the  compensated  imaging  case,  while  we  can  measure  more  differences  in  the  Knox- 
Thompson  (K-T)  case  if  we  want.  This  analogy  is  clear  in  light  of  history,  but  we  can’t  be  too  hard  on 
Knox  and  Thompson  for  not  pointing  it  out.  After  all,  compensated  imaging  did  not  exist  at  the  time 
they  did  their  work. 

Now  in  practice  the  S  (f)  phase  shifts  are  not  perfectly  uniform  over  a  subaperture,  no  matter  how 
small.  Since  we  need  multiple  images  to  get  the  magnitudes  with  Labeyrie’s  method,  we  can  obviously 
use  those  same  images  to  compute  multiple  phase  differences  which  can  be  averaged.  TTie  mean  phase 
difference  will  be  zero,  so  this  process  will  reduce  the  variance.  Since  the  phases  are  b2ised  on  DFTs,  we 
in  fact  need  to  average  them  to  get  reasonable  estimates  for  the  same  reason  Labeyrie’s  method  needs 
multiple  images. 

A  more  serious  problem  occurs  when  the  magnitudes  of  the  spatial  frequencies  become  small  or 
zero.  Then  the  phases  are  either  undefined  or  highly  variable,  and  our  network  of  differences  becomes 
corrupted  or  has  a  hole  rather  than  being  continuous.  TTiis  is  really  the  theoretical  Achilles  heel  of  the 
K-T  method,  and  prevents  it  from  standing  alone  as  a  completely  general  phase  retrieval  scheme.  A 
recent  paper^^  attempts  to  overcome  this  problem. 

We  should  say  a  word  about  the  computations  associated  with  K-T.  If  we  measure  the  phase  dif¬ 
ferences  say  across  every  row  of  our  array  and  just  one  column  (or  vice  versa),  we  have  a  unique  solution 
for  the  phases  which  can  easily  be  computed.  TTiis  was  the  original  approach  of  K-T.  A  rather  obvious 
improvement  was  to  do  it  twice,  once  relying  primarily  on  row  differences,  and  once  relying  on  column 
differences,  and  then  averaging  the  results. 

In  light  of  our  discussion  of  the  compensated  imaging  processing,  we  would  prefer  to  solve  a 
least-squares  problem  using  all  of  the  differences  simultaneously.  Recall  that  the  dimension  of  the 
matrix  problem  matches  the  number  of  subapertures  in  the  compensated  imaging  case.  In  the  K-T  case, 
we  probably  have  10  to  100  spatial  frequencies  in  each  subaperture,  so  rather  than  having  a  least-squares 
problem  of  dimension  say  100,  we  have  one  with  dimension  1000  to  10000  or  more.  Tliis  is  not  cause 
for  panic,  especially  since  we  do  not  need  to  calculate  in  real  time,  but  it  will  require  more  processing. 

Using  the  least-squares  approach  we  could  include  diagonal  differences  as  well,  perhaps  obtaining 
a  more  accurate  result.  Since  the  spatial  frequency  spacing  is  greater  on  the  diagonals  by  V?,  we  might 
assign  those  errors  less  weight.  If  diagonals  are  included,  the  sparse  matrix  being  solved  becomes  less 
sparse.  This  is  not  a  problem,  however,  because  just  as  in  the  compensated  imaging  case  we  would  pre¬ 
invert  the  constant  matrix  A  M  once. 

6.2.2  Triple  Correlation  Method 

If  we  let  the  vector  frequencies  (w,v)  =  f/  and  similarly  the  vector  spatial  coordinates  (jc,y)  =  X/,  we 
can  define  the  triple  correlation  of  an  image  by 


(6.5) 


TTie  Fourier  transform  of  the  triple  correlation  is  known  as  the  bispectrum 

/(3)(f,,f2)  =  /(f,)/(f2)/(-f,-f2). 


(6.6) 
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The  average  bispectrum  of  the  image  can  be  related  to  the  bisp)ectrum  of  the  object  through  the  average 
bispectrum  transfer  function 

,f2)>  =  C><3)(f,  J,)  (5(3)(f|  f^))  (6.7) 

In  our  discussion  of  the  K-T  method  we  have  alluded  to  the  distinction  between  frequency  and 
space  domain  explanations  of  algorithms.  Things  get  much  worse  when  one  contemplates  explaining 
the  triple  correlation  equations  above.  We  have  the  usual  time/space  dichotomy,  but  in  addition  we  can 
group  the  terms  different  ways  to  produce  different  models. 

If  we  group  two  of  the  terms  of  Equation  (6.5)  we  can  consider  them  to  produce  an  autocorrelation 
function  which  is  masked  by  the  third  term.^^*^^  With  certain  types  of  objects  this  can  produce  a  useful 
conceptual  model.  Weigelt^^  adopted  this  speckle  masking  approach  when  he  first  proposed  the  method. 

If  we  consider  the  three  terms  of  Equation  (6.6)  on  an  equal  basis,  we  can  employ  the  phase  clo¬ 
sure  paradigm^^  which  originated  in  radio  astronomy.  The  basic  idea  is  to  add  three  phase  terms  (multi¬ 
ply  three  complex  numbers)  which  have  common  errors  with  alternating  signs  so  that  the  errors  cancel 
out.  Manipulations  of  the  resulting  composite  closure  phases  are  able  to  extract  the  phases  of  interest. 

Naive  application  of  triple  correlation  can  require  prohibitive  computation  and  storage,  so  heuristic 
modifications  of  the  algorithm  have  been  common  with  a  corresponding  lack  of  rigor.  Freeman  et  afi^ 
have  written  an  interesting  paper  comparing  various  triple  correlation  algorithms  to  Knox-Thompson 
with  a  more  appealing  theoretical  basis.  They  also  take  a  frequency  domain  approach  vaguely  similar  to 
our  Knox-Thompson  discussion. 

We  have  defined  the  triple  correlation  and  bispectrum  as  mathematical  objects  without  relating 
them  to  our  problem  of  phase  retrieval.  That  is  accomplished  with  the  aid  of  the  fact^^  that  is  real. 
It  follows  from  Equation  (6.7)  that  the  image  triple  correlation  has  the  same  phase  as  the  object  triple 
correlation. 

We  will  call  the  phase  of  the  bispectrum  P  and  the  phase  of  the  object  ([).  From  Equation  (6.6)  we 
obtain  the  recursive  relationship 

+  f2)  =  4>(fl  )  +  4>(f2)  -  P(fi  J2).  (6.8) 

This  expression  says  that  given  the  phases  at  frequencies  fi  and  (2  and  the  bispectrum  P(f|  f2),  we  can 
compute  the  phase  at  frequency  fi  +  f2.  Recursively  we  can  compute  all  of  the  phases.  To  start  the 
recursion  we  require  a  couple  of  image  phases.  We  can  arbitrarily  set  these  since  they  relate  to  the  posi¬ 
tions^  of  the  object. 

A  key  feature  of  triple  correlation  is  the  multiplicity  of  estimates  of  each  phase  which  are  available. 
We  can  decompose  the  frequency  f  into  many  different  sums  f/+  fy,  each  of  which  yields  a  different  esti¬ 
mate.  While  these  can  be  combined  in  obvious  ways  such  as  averaging,  just  as  in  the  K-T  case  it  is  more 
profitable  to  solve  a  least-squares  problem.  The  minimum-variance  weighting  discussed  by  Freeman  et 
al^^  is  perhaps  the  best  estimator  to  date. 

It  can  be  easily  shown  that  the  four-dimensional  triple  correlation  function  contains  the  K-T  infor¬ 
mation  in  certain  planes.  Triple  correlation  should  always  do  as  well  or  better  than  K-T  since  a  valid 
estimator  could  just  use  the  K-T  data  and  duplicate  the  K-T  algorithm.  By  definition,  a  minimum- 
variance  estimator  which  does  not  use  K-T  must  be  doing  as  well  or  better  than  K-T. 
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The  triple  correlation  algorithms’  biggest  problem  is  their  required  storage  and  processing.  If  we 
have  digital  images  of  dimension  256  by  256,  and  we  store  the  entire  triple  correlation  function  in  single 
precision  4  byte  floating  point,  storage  of  17  Gbytes  is  required,  larger  than  the  largest  magnetic  disks! 
The  processing  required  to  compute  such  a  function  is  also  prohibitive. 

Some  of  the  algorithms  which  estimate  phase  from  the  bisp)ectrum  are  not  too  computationally 
demanding,  but  the  minimum-variance  method  mentioned  above  requires  solution  of  matrices  of  possi¬ 
bly  excessively  large  dimension.  The  demonstration  of  the  method  on  one-dimensional  data  does  not 
give  any  confidence  that  two-dimensional  application  is  feasible. 

This  is  certainly  an  area  of  active  research.  It  is  known^^*^"^  that  there  are  extensive  symmetries  in 
the  triple  correlation.  For  example,  in  the  one-dimensional  case  (two-dimensional  triple  correlation)  we 
can  work  with  one  region  in  the  first  quadrant  which  comprises  about  l/16th  of  the  total  function.  As  the 
theory  progresses,  it  is  hoped  that  the  most  important  regions  of  the  triple  correlation  function  can  be 
identified  and  used  effectively.  The  cited  pap)ers  make  some  progress  in  this  direction. 

6.2.3  Phase  Reconstruction  from  Magnitude 

It  has  been  known  for  years  that  the  problem  of  reconstructing  the  phase  of  a  one-dimensional 
sequence  from  its  magnitude  is  intractable;  there  are  hopelessly  many  solutions.  Until  about  1978 
everyone  assumed  that  the  situation  would  be  no  better  in  the  two-dimensional  case.  The  theoreticians 
did  not  discover  the  error  in  this  assumption  until  practical  algorithms  for  phase  retrieval  were 
published. 

In  1972  Gerchberg  and  Saxton^^  published  an  important  paper  describing  an  iterative  algorithm  for 
phase  retrieval.  Their  application  had  available  magnitude  information  for  both  a  two-dimensional  func¬ 
tion  and  its  Fourier  transform.  For  example,  electron  microscopy  and  X-ray  crystallography  can  satisfy 
this  requirement  by  obtaining  intensities  in  both  the  image  and  diffraction  planes. 

The  algorithm  proceeds  as  shown  in  Figure  6.3.  Initial  random  phases  are  generated,  uniform  on 
the  interval  (-71,71).  They  are  combined  with  the  known  image  magnitudes  and  the  result  is  Fourier 
transformed.  The  phases  of  the  FT  are  likewise  combined  with  the  known  FT  magnitudes  and  the  result 
is  transformed  back.  Again  the  phases  are  combined  with  the  known  image  magnitudes  and  one  itera¬ 
tion  has  been  completed. 

While  interesting,  there  is  no  direct  application  of  this  algorithm  to  phase  retrieval  for  speckle 
imaging.  It  led,  however,  to  the  important  discovery  by  Fienup^^’^^’^^  of  an  iterative  algorithm  suitable 
for  speckle  imaging.  Fienup’s  insight  was  that  one  could  generalize  the  Gerchberg-Saxton  algorithm  so 
that  rather  than  simply  combining  the  FFT  phases  with  the  measured  magnitudes,  one  applied  some 
other  constraint.  Since  we  have  the  FT  magnitudes  in  speckle  imaging,  that  part  of  the  algorithm  can 
proceed  in  the  same  way.  In  the  image  domain,  suitable  constraints  suggested  by  Fienup  include  finite 
support,  in  which  the  image  is  assumed  to  be  zero  outside  some  region,  and  non-negativity,  since  inten¬ 
sities  cannot  be  negative.  Figure  6.4  shows  Fienup’s  algorithm. 

Fienup  supplied  examples,  including  blind  tests,  as  well  as  algorithm  modifications  which  enhance 
convergence.  There  are  certain  ways  the  algorithm  can  stagnate,  and  artifacts  which  may  occur  in  these 
cases.  Fienup  has  given  fixes^^  for  many  of  these  problems.  These  make  interesting  reading  but  they  are 
a  little  beyond  the  scope  of  this  report. 
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Figure  6.3.  Gerchberg-Saxton  algorithm. 


Figure  6.4.  Fienup' s  algorithm. 

Once  Fienup  had  published  his  work,  theorists  went  scurrying  off  to  try  to  explain  why  it  worked. 
The  first  outstanding  issue  concerned  the  uniqueness  or  non-uniqueness  of  the  solutions.  An  algebraic 
explanation  was  finally  produced  involving  polynomial  factoring.  While  any  polynomial  in  one  variable 
can  be  factored,  very  few  polynomials  in  two  variables  have  this  property.  Hence  most  images  have 
unique  phase  solutions.  Hayes^^  has  provided  a  lucid  summary  of  the  mathematics  of  phase 
reconstruction. 
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Once  this  mathematical  foundation  was  laid,  attention  turned  to  better  methods  of  phase  retrieval. 
It  has  been  found  that  while  systems  of  equations  can  be  written  whose  solution  is  the  desired  phase, 
their  solution  is  computationally  impractical  in  all  but  the  smallest  problems. 

The  next  obvious  area  of  mathematical  investigation  is  iterative  algorithms  such  as  Fienup’s  algo¬ 
rithm.  Perhaps  a  more  general  mathematical  understanding  of  such  algorithms  would  lead  to  better  algo¬ 
rithms.  One  mathematical  framework'^®  describes  these  algorithms  in  terms  of  projections  of  convex 
sets  (POCS)  in  Hilbert  space.  An  interesting  recent  paper^*  claims  to  generalize  things  further  by 
appealing  to  some  results  of  nonlinear  programming.  To  make  a  long  story  short,  nothing  has  yet 
emerged  to  challenge  Fienup’s  algorithm  in  the  speckle  imaging  application.  We  do  understand  it  a  lot 
better,  however. 

6.3  SPECKLE  PARAMETERS 

6.3.1  Signal-to-Noise  Ratio 

Speckle  imaging  has  been  developed  and  used  by  astronomers,  so  issues  of  most  concern  to  them 
have  received  the  most  attention.  In  astronomy,  most  of  the  objects  of  greatest  interest  are  extremely 
faint;  bright  objects  are  a  lower  priority.  Faint  objects  also  represent  a  greater  technical  challenge. 

In  space  surveillance,  on  the  other  hand,  in  order  to  obtain  usable  resolution  on  man-made  objects 
using  telescope  ajjertures  on  the  order  of  a  couple  of  meters,  we  consider  mainly  objects  in  relatively  low 
orbits.  For  example,  to  achieve  50  cm  resolution  with  a  2  m  telescope  we  would  confine  our  observa¬ 
tions  to  objects  lower  than  2000  km.  In  this  regime,  with  solar  illumination,  the  satellites  are  quite 
bright  objects. 

It  turns  out  that  the  behavior  of  the  signal-to-noise  ratio  (S/N)  is  different  for  the  two  cases. 
Dainty''^  gives  a  space-domain  calculation  for  the  case  of  a  faint  double  star.  He  computes  the  variance 
of  the  jjeak  correlation,  which  occurs  at  the  binary  star  spacing.  The  main  features  of  his  findings  are  as 
follows: 

(a)  The  rms  S/N  is  proportional  to  the  number  of  detected  photons  per  picture. 

(b)  The  S/N  is  proportional  to  the  telescope  diameter. 

(c)  The  S/N  is  inversely  proportional  to  the  diameter  of  the  long-exjrosure  image,  or 
directly  proportional  to  the  seeing. 

(d)  The  S/N  is  inversely  proportional  to  the  number  of  resolution  cells  in  an  extended 
object. 

Photon  noise  produces  a  bias  in  the  spectral  estimates  which  can  be  removed.  Goodman  and 
Belsher^^  present  a  clear  analysis  of  the  problem  based  on  Poisson  statistics.  It  was  later  shown^^  that  if 
a  photon-counting  camera  is  available,  better  performance  is  obtained  by  making  this  correction  using 
the  actual  number  of  photons  in  the  image  rather  than  the  average  number.  Dainty  and  Greenaway^^ 
also  discuss  noise  reduction  by  elimination  of  frames  containing  less  than  two  photons.  Nisenson  and 
Papaliolios'^  extended  the  analysis  from  simple  speckle  interferometry  to  Knox-Thompson  processing. 
Ayers  et  al  recently  published  an  even  more  comprehensive  papei^^  including  an  analysis  of  triple 
correlation. 
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Miller^^  gives  an  exhaustive  analysis  of  many  noise  sources.  His  assumptions  are  more  general 
than  Dainty’s,  and  his  conclusions  are  somewhat  different  numerically.  It  is  important  to  check  carefully 
just  what  is  being  assumed  in  noise  analyses;  there  are  differences  in  light  levels,  detector  characteristics, 
and  the  quantity  whose  statistics  are  being  specified. 

Let  us  consider  more  carefully  the  case  of  bright  objects.  We  have  already  remarked  above  that 
periodograms  are  rather  poor  estimators.  In  fact,  under  reasonable  (Gaussian)  assumptions  it  can  be 
shown^"^  that  the  standard  deviation  of  periodogram  estimators  is  about  equal  to  the  mean.  This  is  true 
no  matter  how  long  a  data  record  is  Fourier  transformed.  The  explanation  lies  in  the  fact  that  while 
longer  periodograms  involve  more  input  data,  we  are  also  estimating  proportionately  more  parameters. 

It  follows  from  this  that  the  power  spectrum  estimate  of  Labeyrie’s  method  has  S/N  of  approxi¬ 
mately  unity  for  bright  objects.  Note  that  here  we  neglect  detector  noise  along  with  a  few  other  possible 
noise  sources.  These  do  not  change  things  appreciably  in  many  cases  of  imaging  bright  objects  with 
high  quality  equipment.  Hence  unity  is  the  maximum  value  the  S/N  can  assume  in  this  type  of 
processing. 

The  S/N  is  increased  by  averaging  multiple  images.  If  periodograms  of  n  independent  images  are 
averaged,  the  S/N  will  increase  by  a  factor  of  This  can  be  understood  in  many  ways;  for  example, 
the  “signal”  is  being  integrated  coherently  and  its  power  is  increasing  as  while  the  “noise”  is 
integrated  non-coherently  so  its  power  increases  as  n.  The  periodogram  S/N  thus  increases  as  n,  but  we 
take  a  square  root  to  obtain  the  estimate  of  the  spatial  frequency  magnitude  so  its  S/N  increases  as 

If  we  are  interested  in  a  S/N  of  20,  some  4(K)  independent  periodograms  must  be  averaged,  regard¬ 
less  of  the  brightness  of  the  ohjectl  Of  course,  in  the  photon-limited  region  the  individual  S/Ns  are 
lower,  and  sometimes  many  thousands  of  periodograms  are  required. 

This  brief  summary  tried  to  convey  the  idea  that  accurate  S/N  analyses  for  general  cases  are  rather 
complicated.  The  understanding  we  have  achieved  here  will  suffice  for  our  purposes.  A  careful  analysis 
of  S/N  for  each  real  experiment  is  recommended. 

6.3.2  Integration  Time 

So  far  we  have  considered  our  images  as  “instantaneous”  snapshots  through  some  turbulent  atmo¬ 
sphere;  in  reality,  our  detector  must  integrate  for  some  finite  time.  Expressing  this  time  dependence 
Equation  (6.3)  becomes 

/(f,0  =  O(f)5(f,r).  (6.9) 

When  our  detector  integrates  over  a  time  A/,  the  detected  image  is  not  the  /(f)  of  Equation  (6.3)  but 
rather  /(f)  where 

A/ 

/(f)=J/(f,0.  (6.10) 

0 

When  we  do  speckle  processing,  we  have  (ignoring  scaling  factors) 

<l/(f)P>=  l0(f)P<l  j5(f,0P>.  (6.11) 

0 


41 


We  are  interested  in  the  last  term,  the  average  squared  magnitude  of  the  integrated  optical  transfer  func¬ 
tion,  for  various  values  of  At.  Recall  that  Korff'®  studied  this  term  for  long  and  short  exposures. 

Roddier  and  Roddier'^^  generalized  his  results  to  the  intermediate  case  using  the  so-called  “Taylor 
hypothesis.”  When  speckle  patterns  evolve  with  time,  their  behavior  is  usually  considered  to  involve  a 
combination  of  speckle  motion  and  “speckle  boiling,”  in  which  the  individual  speckles  fluctuate.  In  the 
Taylor  hypothesis  the  atmosphere  is  treated  as  a  fixed  snapshot  of  turbulence  which  the  wind  is  blowing 
in  some  direction  at  some  speed.  This  assumption  is  sometimes  questionable,  but  it  did  produce  useful 
results  in  agreement  with  some  experimental  data. 

If  we  consider  the  effect  of  exposure  time  on  the  MTF,  long  exposures  result  in  a  greatly  attenuated 
high-frequency  wing  of  the  MTF,  while  instantaneous  exposures  contain  useful  high-frequency  MTF 
information  out  to  the  diffraction  limit,  as  Korff  showed  theoretically.  In  the  case  of  intermediate  expo¬ 
sure  times,  the  MTF  not  surprisingly  takes  intermediate  values. 

Unfortunately,  Roddier  and  Roddier  required  fourth-order  moments  for  their  calculations,  and  the 
mathematics  is  too  complicated  to  reproduce  here.  We  will  content  ourselves  with  a  few  remarks  on 
their  results. 

For  a  5  m  telescope,  ig  of  13  cm,  and  10  m/s  wind  velocity,  exposure  times  up  to  20  ms  had  only 
negligible  effect  on  the  MTF.  Smaller  telescopes  were  more  affected  by  exposure  time;  under  the  same 
conditions  a  1 .5  m  telescope  suffered  negligible  depression  of  the  MTF  for  exposures  up  to  10  ms. 

High  frequency  components  of  the  MTF  were  more  affected  by  exposure  time  than  low  frequency 
components.  This  would  be  expected  if  we  remember  our  model  of  individual  speckles  moving  around 
the  image.  Each  MTF  frequency  component  can  be  considered  to  be  obtained  by  multiplying  a  spatial 
frequency  basis  function  by  the  image  and  summing.  A  small  motion  affects  the  high  frequencies  more 
than  the  low  frequencies  since  it  involves  motion  of  a  greater  fraction  of  a  spatial  frequency  wavelength 
prior  to  the  correlation  operation. 

In  1978  Karo  and  Schneiderman'^®  published  a  more  comprehensive  and  more  accessible,  although 
perhaps  less  original,  paper  on  exposure  times  and  bandwidths.  They  discuss  the  simulation  of  exposure 
time  by  superimposing  multiple  independent  images.  They  report  measurements  of  the  effect  of  expo¬ 
sure  time  using  the  1 .6  m  telescope  at  the  AMOS  facility  in  Hawaii.  Exposures  of  5  ms  were  the  shortest 
they  used,  but  indications  were  that  shorter  exposures  might  produce  even  better  results.  When  the 
exposure  reached  40  ms,  the  MTF  dropped  into  the  noise  at  about  85  percent  of  the  diffraction  limit. 
Some  MTF  above  the  seeing  limit  was  observed  even  at  320  ms. 

As  long  as  the  MTF  stays  sufficiently  above  the  noise,  correction  for  excessive  exposure  time  is 
theoretically  possible.  Of  course  the  overall  MTF  correction  error  will  increase  as  we  attempt  to  make 
ever  more  heroic  compensations.  This  will  adversely  affect  the  rest  of  our  speckle  imaging  algorithms  in 
some  way.  Hence  it  is  important  to  be  aware  of  the  S/N  when  correcting  for  exposure  time  effects. 

If  one  has  substantially  more  light  per  exposure  than  needed  to  keep  the  single-frame  S/N  near 
unity,  it  clearly  makes  sense  to  shorten  the  exposure  time.  Shortening  the  time  by  a  factor  n  will  increase 
the  total  S/N  by  a  factor  of  as  long  as  the  individual  frames  retain  a  S/N  near  unity. 

On  the  other  hand,  if  the  object  is  very  dim  it  can  be  shown  that  the  single-frame  S/N  is  propor¬ 
tional  to  the  number  of  photons.  Hence  it  pays  to  integrate  longer  and  combine  fewer  frames.  Clearly 
there  is  an  optimum  integration  time  which  is  discussed  by  O’Donnell  and  Dainty.*^^ 
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6.3.3  Bandwidth 


Another  important  issue  is  the  effect  of  optical  bandwidth  on  MTF.  As  Roddier^^  points  out,  there 
are  three  effects  of  optical  bandwidth  which  tend  to  reduce  the  MTF  for  high  spatial  frequencies.  Some 
recent  detailed  theoretical  references^ are  not  the  most  readable  papers  imaginable.  Accounts  of 
good  experimental  work'^^  are  recommended  to  obtain  an  intuitive  grasp  of  bandwidth  effects. 

Perhaps  the  most  basic  effect  of  excessive  bandwidth  is  a  loss  of  temporal  coherence  over  the  tele¬ 
scope  aperture.  Since  the  paths  through  the  turbulence  have  varying  lengths,  causing  the  phase  shift 
which  is  our  main  difficulty,  there  is  a  corresponding  time  delay  for  different  parts  of  the  aperture.  In 
order  to  apply  the  basis  of  speckle  imaging,  the  Van  Cittert-Zemike  theorem,  we  need  the  quasi- 
monochromatic  condition.  We  are  pushing  the  definition  of  quasi-monochromatic  in  practice,  but  we 
certainly  require  some  level  of  temporal  coherence.  Clearly  excessive  bandwidth  can  completely  destroy 
coherence. 

Karo  and  Schneiderman^^  discuss  this  condition,  arriving  at  the  conclusion  that  for  best  results  one 
should  use  Ak/J,  of  0.1  to  0.05.  They  give  a  physical  discussion  using  the  analogy  of  a  good  quality  con¬ 
ventional  lens. 

However  this  should  depend  on  the  telescope  diameter  as  well  as  the  seeing.  It  is  easy  to  see  that 
the  coherence  length  of  a  uniform  spectrum  over  the  bandwidth  AX  is  roughly  7JAX  wavelengths  or 
X  /AX.  (For  example,  if  the  bandwidth  fraction  is  1/10,  the  frequency  extremes  will  be  completely  out  of 
phase  after  about  10  wavelengths.)  The  total  path  length  difference  (phase  shift)  over  the  telescope  aper¬ 
ture  should  not  exceed  this  quantity;  the  problem  is  to  determine  that  phase  shift. 

Assume  that  a  telescope  of  diameter  can  handle  essentially  AX/^  of  60  percent  (entire  visible 
spectrum);  the  wavelength  extremes  are  just  reaching  the  point  of  non-coherence.  If  we  increase  the  tele¬ 
scope  diameter  to  D  >  we  need  to  maintain  coherence  over  the  frequency  extremes  across  the  larger 
diameter  D.  We  can  do  this  by  reducing  the  bandwidth  proportionately,  arriving  at  the  often-quoted  rule 


The  second  bandwidth  effect  is  the  dispersion  introduced  by  the  atmosphere  if  the  telescope  is  not 
pointing  directly  at  the  zenith.  This  can  be  corrected  by  a  prism  arrangement  and  so  represents  less  of  a 
fundamental  limitation  than  the  other  effects. 

The  third  effect  is  a  wavelength-dependent  magnification  of  the  entire  image.  As  discussed  by 
Roddier,^^  this  leads  to  a  requirement  similar  to  Equation  (6.12). 

The  situation  is  quite  similar  to  the  exposure  time  case.  It  is  necessary  to  push  the  bandwidth  as 
high  as  possible  to  enhance  S/N  ratio.  When  this  is  done,  the  higher  portion  of  the  MTF  is  depressed. 
Again,  it  is  theoretically  possible  to  correct  for  excessive  bandwidth  subject  to  ultimate  S/N  limitations. 
Just  as  in  the  exposure  time  case,  the  theory  involves  enough  approximations  to  keep  one  alert  to  the 
possibility  that  experimental  conditions  may  cause  unexpected  behavior.  Any  available  experimental 
measurements  of  bandwidth  effects  during  data  collection,  or  with  the  equipment  used  for  data  collec¬ 
tion,  are  probably  worth  examining. 
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6.3.4  Field  of  View 


The  field  of  view  over  which  speckle  imaging  or  adaptive  optics  works  is  often  loosely  defined  to 
be  the  “isoplanatic  patch”  which  variously  is  described  as  1”  to  5”  or  even  more.  Fried  wrote  an  often 
referenced  paper^"^  providing  theoretical  analysis  of  five  different  types  of  isoplanicity,  one  of  which  was 
speckle  interferometry  isoplanicity.  The  differences  essentially  involve  the  equipment  used  and  the 
exposure  times.  The  goal  was  to  characterize  isoplanicity  in  each  case  so  that  corrections  could  be  made 
for  its  effects.  Another  important  theoretical  paper^^  by  Korff  et  aL  gives  a  deeper  theoretical  analysis 
of  isoplanicity  applicable  for  the  speckle  imaging  case.  This  is  not  the  easiest  paper  in  the  world  to  read. 

It  is  again  more  helpful  for  our  purposes  to  look  at  some  exj>erimental  work.  Schneiderman  and 
Karo  discussed  the  results  of  some  measurements  at  AMOS  in  1976  in  a  paper^^  on  binary  star  speckle 
interferometry.  Unfortunately,  they  specialized  their  analysis  to  the  problem  of  resolving  two  (individu¬ 
ally)  unresolvable  stars.  The  resulting  parametric  approach  is  based  on  fringe  visibility  and  is  not 
directly  applicable  to  extended  objects.  They  also  had  some  experimental  difficulties  which  cast  some 
doubt  on  their  experimental  results.  The  importance  of  their  paper  is  that  it  used  the  technique  of  com¬ 
puting  the  cross-correlation  between  two  speckle  images  formed  from  double  stars  of  known  separation. 
Another  similar  paper^^  on  measures  of  isoplanicity  using  binary  star  systems  was  published  at  about 
the  same  time. 

The  cross-correlation  gives  us  a  measure  of  the  nonisoplanatic  effect.  If  we  take  out  the  average 
separation  of  the  two  images  a  and  /?,  and  let  the  coordinate  system  of  each  image  be  centered,  the 
simple  cross-correlation  mentioned  above  becomes  one  point  of  the  cross-correlation  function  If 

we  consider  the  more  general  correlation  function  RahOO  where  x  is  the  magnitude  of  vector  correlation 
distance,  we  study  its  Fourier  transform  the  cross-spectrum  The  integral  of  S^b  is  /?j/,(0).  Korff 

et  showed  the  important  role  of  the  cross-spectrum  as  a  measure  of  the  effect  of  nonisoplanicity  on 
speckle  imaging.  As  might  be  expected,  the  higher  spatial  frequencies  are  attenuated  more  in  the  cross¬ 
spectrum. 

Ebersberger  and  Weigelt  published  an  interesting  study^^  in  which  they  computed  correlation 
functions  and  cross-spectra  of  speckle  images  of  binary  stars  with  various  separations.  The  correlation 
ranged  from  55  percent  for  separations  of  2.9”  down  to  8  percent  for  7.8”  separations.  In  addition  they 
showed  cross-spectra  of  binary  stars. 

For  2.9”  separation,  the  cross-spectrum  dropped  to  0.5  at  a  normalized  spatial  frequency  u  of 
about  0.6;  it  dropped  to  zero  at  about  0.9.  For  the  4.3”  separation  the  drop  to  0.5  occurred  at  frequency 
0.35,  with  zero  again  reached  at  0.9.  For  the  7.8”  separation  the  numbers  were  about  0.1  and  0.8. 

Christou  and  Hege^^  continued  in  this  vein,  showing  additional  plots  of  normalized  cross¬ 
spectrum.  They  identified  this  function  as  the  correction  factor  required  in  speckle  holography. 

We  have  not  discussed  this  technique,  but  briefly  if  a  point  source  is  near  the  object  being  imaged, 
the  phase  of  the  object  may  be  retrieved  directly.  From  an  autocorrelation  point  of  view^  if  the  point 
source  is  at  a  distance  of  more  than  twice  the  diameter  of  the  object,  the  autocorrelation  function  of  the 
image  will  consist  of  three  parts.  The  part  near  the  origin  may  be  discarded  and  the  cross-terms  will  be 
images  of  the  object.  This  technique  also  can  be  explained^^  by  analogy  to  holography. 
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Reluming  to  our  discussion  of  isoplanicity,  since  the  reference  object  in  speckle  holography  has 
some  angular  separation  from  the  object  being  imaged,  we  have  a  need  for  a  correction.  Dividing  by  the 
normalized  cross-spectrum  effects  this  correction. 

We  can  interpret  drops  in  the  cross-spectrum  as  reductions  in  S/N.  Hence  we  lose  the  most  S/N  at 
the  high  spectral  frequencies.  In  speckle  imaging,  if  we  have  enough  data  of  good  S/N,  we  can  over¬ 
come  this  problem  by  statistical  averaging.  If  we  don’t,  we  will  lose  some  resolution.  Things  seem 
worse  for  adaptive  optics.  The  instantaneous  estimate  of  the  turbulence  phase  shift  will  be  noisier,  so 
our  real-time  correction  will  be  poorer,  especially  at  the  high  spatial  frequencies. 

In  conclusion,  isoplanicity  is  an  elusive  concept  which  must  be  considered  in  the  context  of  the 
actual  system  under  study  rather  than  as  an  absolute  quantity  def)endent  only  on  the  atmosphere.  It 
seems  inescapable  that  its  effects  will  be  observed  whenever  the  field  of  view  of  an  image  exceeds  one  or 
two  arc  seconds.  It  may  represent  a  substantial  problem  for  satellite  imaging  due  to  the  large  angular 
extent  of  many  low-eanh-orbit  satellites. 
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7.  APPLICATION  TO  SPACE  SURVEILLANCE 


7.1  INTRODUCTION 

We  have  completed  our  review  of  the  new  adaptive  and  speckle  techniques  for  high-resolution 
optical  imaging.  By  this  time,  the  reader  should  be  able  to  follow  the  gist  of  the  literature  on  these  sub¬ 
jects.  To  that  end,  we  have  used  much  of  the  terminology  and  notation  which  is  current  with  researchers 
in  these  areas. 

We  wish  to  consider  the  new  techniques  together  in  this  section,  taking  a  more  unifying  viewpx)int. 
The  adaptive  system  will  be  analyzed  as  a  coherent  integration  prior  to  the  noncoherent  speckle 
processing. 

We  will  also  discuss  the  peculiar  problems  of  imaging  satellites  rather  than  the  astronomical 
objects  which  dominate  the  unclassified  literature.  We  will  consider  the  limits  on  integration  time 
imposed  by  satellite  motion.  The  bistatic  nature  of  incoherent  optical  imaging  will  be  contrasted  with 
typically  monostatic  alternatives  such  as  radar  imaging. 

Adaptive  systems  will  always  contain  small  errors.  Since  we  are  adaptively  correcting  phase,  these 
errors  will  manifest  themselves  in  this  domain.  Atmospheric  scintillation  will  also  produce  small  errors 
which  are  unaffected  by  adaptive  phase  corrections.  We  will  explain  the  similar  effects  of  small  ampli¬ 
tude  and  phase  errors  on  the  OTF  by  appealing  to  well-known  results  of  communications  theory. 

To  make  the  transition  from  speckle  interferometry  to  speckle  imaging  we  need  to  select  a  phase- 
retrieval  algorithm.  We  will  provide  guidelines  for  that  selection  decision  which  depend  on  the  data;  our 
satellite  imaging  application  will  then  dictate  the  decision. 

We  have  discussed  the  various  errors  and  degradations  which  can  hinder  the  new  imaging  tech¬ 
niques  at  some  length.  While  it  is  obviously  preferable  to  avoid  these  errors  and  degradations  as  much 
as  possible  during  experiment  design  and  data  collection,  there  are  times  when  unknown  degradations 
have  affected  our  data.  Most  astronomical  data  collections  may  be  repeated  at  our  leisure  (excepting  the 
occasional  supernova)  subject  only  to  equipment  availability.  This  is  not  the  case  in  satellite  surveil¬ 
lance,  since  satellite  lifetimes  are  comparatively  short  and  the  most  interesting  data  may  involve  short, 
unique  events. 

In  cases  where  a  data  set  is  irreplaceable,  it  is  desirable  to  develop  techniques  of  evaluating  its 
quality  as  well  as  its  potential  for  enhancement  using  speckle  techniques.  The  interferometer  signals 
used  in  adaptive  optics  are  one  clear  source  of  data  quality  measurement,  but  they  do  not  exist  in  the 
speckle  case  and  may  not  always  be  recorded  in  adaptive  systems  due  to  the  volume  of  data.  We  will 
describe  some  measures  of  data  quality  which  are  calculated  from  short-term  images  themselves. 

7.2  GENERALIZED  OPTICAL  SURVEILLANCE 
7.2.1  Block  Diagram 

For  purposes  of  discussion,  consider  the  generalized  system  depicted  in  Figure  7.1.  The  dashed 
boxes  may  or  may  not  exist  in  a  given  system.  This  may  be  thought  of  either  as  the  combination  of 
adaptive  optics  and  optional  speckle  post-processing,  or  speckle  imaging  with  optional  adaptive 
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pre-processor.  The  object  being  imaged  is  a  satellite.  We  allow  the  pMjssibility  of  using  a  reference 
object  to  assist  either  the  adaptive  optics  or  the  speckle  processing.  This  might  be  either  a  natural  refer¬ 
ence  object  such  as  a  star  or  artificial  stars  created  by  a  laser^^  in  the  atmosphere. 
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Figure  7.7 .  Generalized  optical  sun'eillance  system. 


There  are  two  physical  effects  that  we  are  estimating  no  matter  how  we  operate  this  system.  One  is 
some  measure  of  turbulence,  either  on  a  short-term  or  long-term  basis  or  both.  The  other  is  of  course  an 
estimate  or  sequence  of  estimates  of  the  spatial  frequencies  of  the  target  itself.  The  former  is  used  to  cor¬ 
rect  the  latter. 

7.2.2  Effect  of  the  Wavefront  Modifier 

Light  from  the  satellite  being  imaged  passes  through  the  turbulent  atmosphere  and  enters  the 
optional  wavefront  modifier.  This  device  can  be  considered  to  be  a  time-varying  matched  filter.  It  pro¬ 
vides  coherent  integration  gain  by  matching  the  effect  of  the  turbulence. 

Recall  from  communications  theory  that  when  a  signal  with  some  specified  frequency  distribution 
is  received  in  the  presence  of  white  noise,  a  matched  filter  maximizes  the  output  S/N.  The  matched  filter 
has  frequency  spectrum  which  is  the  complex  conjugate  of  the  signal  sp>ectrum,  or  equivalently  impulse 
response  which  is  the  time-reversed  signal  itself. 
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There  are  some  differences  in  the  physics  of  our  situation  but  very  little  difference  in  the  mathe¬ 
matics.  Rather  than  time  functions  we  are  concerned  with  space  functions.  The  finite  waveform  dura¬ 
tion  T  of  matched  filter  derivations  is  replaced  by  the  telescope  aperture  D.  We  also  are  working  in  two 
spatial  dimensions  rather  than  one  time  dimension;  the  two-dimensional  analysis  is  a  straightforward 
extension  of  the  one-dimensional. 

Since  we  are  performing  our  wavefront  correction  in  the  pupil  plane,  we  are  dealing  with  the 
mutual  coherence  function  so  it  makes  sense  to  think  of  the  frequency-domain  interpretation  of  the 
matched  filter.  The  atmosphere  is  being  modeled  as  a  filter  with  unity  magnitude  at  each  spatial  fre¬ 
quency  and  some  phase  shift.  Each  subaperture  corrects  for  the  atmospheric  phase  shift  by  producing  a 
conjugate  phase  shift,  the  definition  of  the  matched  filter. 

What  remains  is  to  quantify  the  coherent  integration  gain  that  we  obtain  from  a  perfect  wavefront 
matched  filter.  Recall  the  frequency  domain  heuristic  explanation  of  speckle  imaging  discussed  previ¬ 
ously  in  Section  6. 1.1.1.  Based  on  the  interferometric  view  of  imaging,  we  can  imagine  the  transfer 
function  of  the  atmosphere  to  be  composed  of  a  two-dimensional  array  of  subapertures  whose  phase  shift 
is  constant  and  independent  of  the  other  subapertures.  We  can  use  this  model  both  to  compute  the 
transfer  function  5(f)  as  before  and  to  derive  the  S/N  advantage  of  the  matched  filter. 

We  previously  showed  that  the  speckle  transfer  function  grows  as  the  pupil  redundancy  A^(f)  rather 
than  as  N^{f)  due  to  the  noncoherent  nature  of  Labeyrie  interferometry.  This  implies  that  the  spatial  fre¬ 
quency  MTF  is  proportional  to  VN(f). 

Unfortunately,  we  cannot  directly  achieve  the  theoretical  S/N  of  the  speckle  process  but  must 
instead  compute  periodograms  which  have  poor  statistical  properties,  e.  g.,  a  standard  deviation  on  the 
order  of  the  mean.  Hence  we  face  a  S/N  limit  at  all  spatial  frequencies  of  unity.  This  could  represent  a 
tremendous  loss  in  the  case  of  very  high  light  levels  and  fixed  integration  time.  Practically  we  have 
finite  light,  and  we  are  able  to  select  our  integration  time  if  we  use  suitable  equipment.  The  solution  is  to 
use  very  short  integration  intervals  and  many  of  them,  so  that  the  photon-limited  S/N  per  frame  is  not 
much  higher  than  unity  to  begin  with. 

In  the  case  of  a  matched  filter  wavefront  corrector,  the  MTF  is  just  the  magnitude  of  the  perfect  cir¬ 
cular  aperture  transfer  function  which  is  proportional  to  N(f). 

Now  that  we  understand  qualitatively  the  effect  of  coherent  and  noncoherent  integration  on  the 
MTF,  we  will  turn  our  attention  to  the  S/N  improvement  due  to  matched  filtering.  If  we  consider  the 
image  as  a  whole,  we  will  obtain  a  coherent  integration  gain  of  D/r^  in  each  direction,  or  a  total  integra¬ 
tion  gain  of  (Dlr^)^  over  a  telescope  of  aperture  r^. 

We  are  more  interested,  however,  in  the  gain  at  individual  spatial  frequencies  fjt.  Recall  that  for 
bright  images  the  maximum  S/N  of  a  periodogram  frequency  is  unity  for  frequencies  beyond  the  seeing 
limit.  The  matched  filter  circumvents  this  limitation  by  coherently  integrating  spatial  frequencies  all  the 
way  to  the  diffraction  limit.  Hence  the  seeing  limit  becomes  the  diffraction  limit  and  the  S/N  can  exceed 
unity  at  all  spatial  frequencies.  Of  course  the  improvement  at  a  given  frequency  f  will  depend  on  N{f), 

Numerically  we  are  motivated  to  find  the  transfer  function  7(0  of  an  ideal  circular  aperture  which 
Equation  (6.2)  relates  to  N{f).  We  have  defined  f  =  (w,v)  but  T{f)  actually  only  depends  on  the  magni¬ 
tude  of  this  vector  frequency,  .  We  are  interested  in  computing  the  ratio  of  the  overlap  area 

circles  in  Figure  3.5  to  the  area  of  one  circle.  A  straightforward  calculation  provides  us  with  the  result 
that 
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where  g  is  the  ratio  of  the  vector  spatial  frequency  magnitude  to  the  spatial  frequency  diffraction  limit. 
This  function  is  shown  in  Figure  7.2. 


Figure  7.2.  Circular  aperture  transfer  function. 


We  can  obtain  the  effect  of  this  matched  filter  from  the  graph  for  a  given  scenario.  For  example, 
assume  that  we  have  a  1.0  m  telescope,  and  that  rg  is  10  cm.  Then  A'(O)  is  {Dlrf)^  or  IfX).  Figure  7.3 
depicts  some  idealized  curves  for  various  situations.  In  the  lower  left,  we  start  with  a  10  cm  telescope 
imaging  some  object  with  maximum  S/N  of  some  arbitrary  value  “R.”  If  we  increase  our  aperture  to 
1  m,  we  obtain  a  S/N  which  is  100^  (40  dB)  greater  for  the  0  spatial  frequency.  These  two  curves  are 
just  rescaled  logarithms  of  Figure  7.2. 

The  top  curve  in  Figure  7.3  illustrates  perfect  integration  and  is  just  a  stretched  version  of  the 
simple  1  m  telescope  curve.  The  theoretical  speckle  curve  shows  that  below  the  seeing  limit  coherent 
combination  occurs  and  we  follow  the  simple  1  m  telescope  curve.  The  spatial  frequencies  beyond  the 
seeing  limit  (i.  e.,  normalized  frequency  0.1)  are  noncoherently  integrated.  Note  the  resemblance  to  the 
speckle  MTFs  obtained  by  Korff  and  depicted  in  Figure  6.2.  Unfortunately  we  must  estimate  these  spa¬ 
tial  frequencies  using  periodograms,  which  reduce  the  S/N  to  unity  for  a  single  image  as  shown  in  the 
bottom  curve.  When  the  seeing  limit  is  exceeded,  higher  frequencies  are  limited  to  unity  S/N  until  we 
are  so  near  the  diffraction  limit  that  further  decreases  occur. 

The  theoretical  speckle  curve  is  obtained  by  the  following  construction.  Since  the  integration  gain 
of  the  noncoherent  speckle  process  is  the  square  root  of  the  coherent  gain  of  the  1  m  perfect  integration 
curve,  we  take  half  this  distance  on  the  log  scale.  Since  we  have  an  offset  of  R  dB,  the  speckle  curve  lies 
halfway  between  the  top  curve  and  the  line  y  =  R.  Where  the  simple  1  m  telescope  curve  exceeds  this 
value,  we  follow  that  curve  instead. 
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Figure  7.3.  Effects  of  coherent  integration. 

Let  us  take  a  numerical  example,  say  at  g  of  0.5,  well  beyond  the  seeing  limit.  Tig)  here  is  0.39 
which  implies  that  N(g)  is  3900  (36  dB)  in  the  coherent  case  and  the  square  root  of  this  quantity  (62  or 
18  dB)  in  the  noncoherent  case.  The  advantage  of  coherent  integration  is  evident. 

The  loss  in  S/N  in  going  from  the  theoretical  speckle  curve  to  the  periodogram  speckle  case 
depends  on  the  value  of  R.  If  R  is  small  enough  (below  0  dB)  there  may  be  very  little  loss.  This  is  syn¬ 
onymous  with  short  integration  times,  which  are  desirable  for  speckle  anyway  since  they  provide  the 
best  speckle  transfer  function. 

7.2.3  Adaptive  Wavefront  Sensing 

If  an  adaptive  matched  filter  is  being  used,  obviously  real-time  wavefront  sensing  and  processing 
are  needed.  The  wavefront  detector  shown  in  Figure  7. 1  would  presumably  split  off  some  light  to  an 
interferometer  arrangement  which  would  then  be  converted  to  wavefront  modifier  drive  signals  in  the 
processor.  Section  5  has  adequately  discussed  compensated  imaging.  We  will  make  a  few  remarks 
about  new  directions  in  adaptive  systems  and  their  possible  application. 

Compensated  imaging  typically  uses  light  reflected  from  the  target  for  turbulence  sensing.  An 
alternative  theoretically  permitting  the  imaging  of  fainter  objects  is  to  use  some  bright  reference  source. 
A  bright  star  is  the  first  thing  that  comes  to  mind,  and  indeed  one  is  shown  in  Figure  7.1,  but  there  is  an 
extreme  shortage  of  such  stars. 
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An  alternative  is  to  use  a  so-called  artificial  star  created  by  directing  laser  light  at  certain  layers  of 
the  atmosphere.  For  example,  experiments  have  been  performed^^  using  the  sodium  layer  at  about 
90  km  altitude.  Adaptive  optics  require  a  reference  which  is  unresolvable  by  a  subaperture.  Since  a 
subaperture  is  on  the  order  of  in  diameter,  Thompson  and  Gardner^^  maintain  that  the  angular  extent 
of  an  artificial  star  for  this  purpose  should  be  about  1  arc  sec.  Smaller  extents  produce  less  light,  degrad¬ 
ing  adaptive  performance  due  to  simple  S/N  considerations.  Larger  extents  degrade  adaptive  perfor¬ 
mance  as  the  reference  is  partially  resolved  by  the  subapertures. 

Future  adaptive  optics  systems  will  likely  operate  digitally  from  the  interferometer  intensity  detec¬ 
tors  through  the  processor.  D/A  converters  will  then  drive  the  wavefront  modifier  circuitry.  Digital 
operation  implies  segmentation  of  time  rather  than  continuous  calculation  as  in  the  analog  approach. 
This  segmentation  is  compatible  with  an  operating  mode  in  which  a  reference  star  is  switched  on  and  off 
and  the  wavefront  sector  synchronously  samples  data  only  when  reference  light  is  entering  the  instru¬ 
ment.  Such  a  mode  will  probably  be  necessary  to  avoid  contamination  of  the  image  data  by  the  broad¬ 
band  reference  star  light. 

7.2.4  Spectral  Filter 

The  adaptive  optics  components  depicted  in  Figure  7.1  are  assumed  to  work  in  white  light.  Hence 
it  seems  logical  to  image  in  white  light  as  well.  However,  we  are  assuming  that  speckle  techniques  may 
be  used  to  post-process  the  data  and  enhance  it.  We  thus  need  to  satisfy  the  speckle  bandwidth  restric¬ 
tions  discussed  in  Section  6.3.3. 

The  exact  bandwidth  of  this  filter  is  problematical.  We  have  suitable  guidelines  for  non-adaptive 
systems,  although  the  optimum  bandwidth  does  depend  on  the  seeing.  In  this  application  optimum 
bandwidth  depends  on  the  seeing  and  also  on  the  performance  of  the  adaptive  components  of  the  instru¬ 
ment.  The  adaptive  equipment  in  some  cases  may  be  considered  to  convert  a  small  fg  into  a  much  larger 
one,  permitting  a  substantial  relaxation  of  the  bandwidth  restrictions. 

The  ideal  situation  would  be  to  divide  the  spectrum  into  multiple  bands,  recording  each  band.  One 
could  then  make  the  decision  during  the  data  reduction  process.  Experience  with  a  real  system  would 
doubtless  lead  to  guidelines  permitting  suitable  filter  selection  even  without  multiple  spectral  bands. 

7.2.5  Image  Detector 

New  photon-counting  detectors^^’ are  becoming  available.  These  have  major  advantages  for  use 
in  collection  of  data  which  may  be  enhanced  with  speckle  techniques.  The  most  obvious  advantage  is 
that  the  integration  interval  can  be  determined  during  data  processing  rather  than  before  collection.  This 
gives  the  user  the  chance  to  optimize  this  parameter  in  every  case;  it  is  particularly  useful  in  the  satellite 
imaging  case  where  we  have  already  noted  that  unique  events  are  common. 

Another  obvious  advantage  of  photon-counting  detectors  is  that  the  data  are  recorded  digitally. 
This  saves  a  lot  of  tedious  and  expensive  digitization.  For  dim  objects  the  volume  of  data  can  also  be 
reduced. 

We  have  noted  that  it  is  advantageous  to  select  bandwidth  and  integration  time  after  the  data  are 
collected.  It  seems  that  what  is  really  needed  is  a  color  photon-counting  camera  which  records  photon  x 
and  y  coordinates,  arrival  times,  and  some  measure  of  wavelength.  This  would  provide  maximum 
flexibility  and  the  potential  for  either  true  color  images  or  improved  S/N. 
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The  method  of  producing  color  images  is  obvious.  Improving  S/N  relies  on  the  likelihood  that 
most  objects  have  similar  spatial  frequency  distribution  over  a  substantial  range  of  wavelengths.  If  this 
is  the  case,  the  S/N  could  be  increased  by  combining  images  from  various  color  bands.  A  scaling  correc¬ 
tion  for  the  magnification  factor  would  probably  be  necessary. 

7,2.6  Post-Processing 

Any  system  has  its  imperfections,  and  despite  our  analysis  of  the  wavefront  corrector  as  a  perfect 
matched  filter,  there  will  be  certain  degradations  in  the  image  received  by  the  detector.  For  example,  the 
photon  noise  associated  with  the  finite  number  of  photons  received  in  each  subaperture  will  produce 
some  error  in  the  differential  wavefront  tilts.  The  calculation  of  the  wavefront  modifier  weights  will 
involve  some  error. 

This  adaptive  system  is  a  feedback  control  system.  All  such  systems  work  by  sensing  errors  at  the 
output  of  the  system  and  driving  the  system  in  such  a  way  as  to  reduce  those  errors.  The  magnitude  of 
the  residual  error  depends  on  the  feedback  loop  gain,  but  the  point  is  that  zero  error  is  never  achieved. 

There  are  also  physical  effects  such  as  scintillation.  We  are  modifying  only  the  phase  of  the 
received  wavefront,  so  errors  in  the  amplitude  remain,  although  those  are  of  course  much  smaller.  We 
discuss  these  errors  in  more  detail  below, 

A  more  fundamental  problem  is  the  problem  of  isoplanicity.  An  adaptive  system  containing  a 
single  wavefront  modifier  operates  correctly  only  over  some  isoplanatic  angle  of  perhaps  several  arc  sec¬ 
onds.  Targets  which  subtend  more  than  this  angle  will  be  imaged  with  errors  contributing  to  a  loss  of 
resolution. 

Much  of  the  theory  of  speckle  imaging  looks  applicable  to  the  problem  of  post-processing  data 
which  has  been  adaptively  processed.  There  will  be  differences  between  straight  speckle  processing  and 
speckle  post-processing  of  adaptively  corrected  data,  and  theoretical  work  remains  to  be  done. 

If  adaptive  errors  can  be  statistically  characterized,  combining  multiple  images  in  the  speckle  tradi¬ 
tion  seems  potentially  able  to  reduce  those  errors  and  produce  improved  images.  This  approach  seems 
more  fruitful  and  more  true  to  the  physics  than  the  blind  application  of  image  processing  techniques  to 
filter  out  noise  or  enhance  edges.  Wiener  filtering  of  individual  frames  is  likely  to  be  inferior  to  tech¬ 
niques  based  on  multiple  frames.  In  the  case  of  a  lack  of  isoplanicity,  Wiener  filtering  is  not  even  theo¬ 
retically  applicable. 

The  isoplanicity  problem  is  one  area  in  which  post-processing  may  improve  things.  There  is  no 
obvious  way  to  deal  with  this  problem  in  the  adaptive  framework  with  a  single  wavefront  modifier  and 
still  image  objects  subtending  large  angles.  However,  one  might  be  able  to  divide  a  speckle  image  into 
sectors,  performing  speckle  processing  on  each  sector  individually. 

In  sum,  there  is  no  shortage  of  possibilities  for  post-processing.  We  will  discuss  some  of  these 
ideas  further  later  in  this  section. 

73  SPECIAL  PROBLEMS  OF  SATELLITE  IMAGING 

We  have  already  alluded  to  some  of  the  problems  of  satellite  imaging  as  compared  to  astronomical 
imaging.  Next  we  will  discuss  these  in  more  detail. 
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7.3.1  Satellite  Motion 


Low-earth-orbit  satellites  suitable  for  optical  imaging  typically  pass  over  the  telescope  rather  rap¬ 
idly.  One  implication  of  this  is  that  high  performance  trackers  are  needed  to  keep  the  telescope  pointed 
at  the  target.  We  will  not  consider  that  issue  any  further,  but  it  is  a  major  hardware  consideration. 

Another  implication  which  concerns  us  is  the  limited  length  of  the  available  image  sequence. 
Regardless  of  the  satellite  motion,  during  the  twilight  hours  suitable  for  imaging  we  expect  to  collect 
data  for  only  several  minutes.  This  is  substantially  less  than  the  tens  of  minutes  or  hours  often  used  in 
speckle  data  collection.  We  have  a  potential  S/N  problem;  the  coherent  integration  provided  by  wave- 
front  compensation  could  be  mandatory  here  for  imaging  certain  objects. 

Modem  high-resolution  optical  techniques  are  quite  effective  in  astronomy  in  which  most  objects 
of  interest  appear  to  be  two-dimensional  objects  attached  to  the  “celestial  sphere.”  Rotations  about  the 
line-of-sight  (LOS)  are  easily  removed  from  the  data. 

In  satellite  imaging,  we  image  three-dimensional  objects  whose  rotation  in  general  is  not  about  the 
LOS.  We  cannot  correct  for  the  component  of  motion  not  about  the  LOS  so  unavoidable  blurring  is  pro¬ 
duced.  This  effect  further  limits  the  length  of  data  which  we  can  integrate. 

We  would  like  to  know  how  long  we  can  expect  to  integrate  data  for  speckle  processing.  We  are 
trying  to  gain  resolution  through  speckle  integration  but  we  are  losing  some  resolution  through  blurring. 
We  can  calculate  this  effect  for  a  given  example  with  some  simplifying  assumptions  to  provide  insight 
into  the  problem. 


Imaging  Interval 


Consider  the  geometry  of  Figure  7.4.  Assume  that  we  are  imaging  a  satellite  in  a  circular  orbit  of 
3(X)  km  which  passes  directly  over  the  observing  site.  The  orbital  period  is  90.5  min,  which  gives  a  tar¬ 
get  velocity  of  about  7.7  km/s.  When  the  target  is  directly  overhead,  the  angular  rate  of  change  is 
1.487s. 

We  will  image  starting  at  the  zenith  for  some  interval  as  shown.  The  satellite  being  imaged  con¬ 
sists  of  a  featureless  cylinder,  the  cross-section  of  which  is  shown  in  the  figure,  and  a  square  solar  panel 
inclined  at  some  angle  0  with  respect  to  broadside.  We  are  interested  in  high-resolution  imaging  of  this 
solar  panel.  Even  if  the  panel  is  tracking  the  sun,  we  expect  very  little  angle  change  over  the  short  imag¬ 
ing  interval,  so  we  will  assume  that  0  changes  only  due  to  the  change  in  viewing  angle  from  the  ground. 

Assume  further  that  we  are  using  a  2  m  telescope,  which  at  SOO-nm  wavelength  gives  us  a 
diffraction-limited  resolution  of  about  0.05  second  of  arc.  Let  the  field  of  view  of  the  telescope  be  10  m. 
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and  the  dimension  of  the  solar  panel  be  5  m.  Assume  that  the  blurring  effect  is  due  to  foreshortening  of 
the  panel  as  its  angle  changes  in  the  x  direction,  as  shown  in  Figure  7.5.  Note  that  there  is  no  blurring  in 
the  y  direction  in  this  example.  We  assume  in  Figure  7.5  that  0  is  40“  at  zenith,  and  increases  to  50*  over 
a  potential  10“  imaging  interval. 


(a)  Panel  at  zenith  (0  =  40“) 


Figure  7.5.  Imaging  a  solar  panel. 


The  effect  of  this  foreshortening  on  imaging  is  easy  to  determine  approximately.  If  the  effective 
width  of  the  panel  changes  10  percent  over  the  imaging  interval,  then  spatial  frequencies  beyond  about 
10  cycles  over  the  panel  are  lost.  At  that  point,  things  are  completely  out  of  phase  from  one  side  of  the 
panel  to  the  other. 

In  our  example,  the  angular  extent  of  the  solar  panel  at  zenith  is  about  (5  m)(cos  40“)/(300  km)  or 
2.6”.  The  diffraction-limited  maximum  spatial  frequency  is  then  about  2.6/.05  or  52  cycles/panel. 
Clearly  we  can  tolerate  something  like  2  percent  foreshortening. 

Figure  7.6  shows  curves  of  integration  time  vs  spatial  frequency  limit  measured  as  a  fraction  of  the 
diffraction  limit.  Remember  that  the  diffraction  limit  in  cycles/meter  for  the  panel  dep)ends  on  its  orien¬ 
tation,  being  a  maximum  at  broadside.  We  can  see  that  if  we  confine  our  integration  to  about  0.85  s  blur¬ 
ring  will  not  limit  our  resolution.  If  we  go  to  2  s  to  improve  the  S/N,  we  lose  more  than  half  our  poten¬ 
tial  resolution.  The  situation  gets  better  for  smaller  0  as  shown;  it  also  gets  worse  for  larger  0. 

An  actual  satellite  is  more  complicated  than  the  simple  model  analyzed  here.  It  is  likely  in  practice 
that  many  different  surfaces  of  the  satellite  will  be  at  different  viewing  angles.  Hence  blurring  will 
affect  different  surfaces  differently.  This  type  of  analysis  can  be  used  in  interpreting  the  results  of 
integrated  optical  imagery  by  considering  portions  of  the  object.  This  is  particularly  easy  if  a  good 
motion  solution  is  available. 

13.2  Bistatic  Imaging 

In  the  case  of  imaging  distant  astronomical  objects,  we  may  consider  them  fixed  to  the  celestial 
sphere.  That  is,  we  are  at  the  center  of  the  celestial  sphere  so  a  given  object  always  looks  the  same. 


55 


10 


r(s)  3 


5 


2 


0 


0.2 


0.4 


0.6 


0.8 


1.0 


Normalized  Spatial  Frequency  Magnitude 


Figure  7.6.  Maximum  integration  time  for  specified  resolution. 


except  for  a  rotation.  Hence  a  single  parameter  describes  the  family  of  images  that  we  can  obtain.  We 
also  know  what  the  rotation  is,  so  it  is  easy  to  correct. 

If  the  receiver  and  transmitter  of  an  imaging  system  are  co-located,  we  are  performing  monostatic 
imaging.  This  is  the  case  in  radar,  for  example.  When  imaging  a  three-dimensional  object  we  can  char¬ 
acterize  the  possible  views  of  that  object  by  imagining  the  object  to  be  at  the  center  of  a  large  sphere.  As 
we  move  around  the  outside  of  the  sphere,  we  obtain  different  projections  of  the  object.  Hence  two  addi¬ 
tional  parameters  “latitude  and  longitude’’  join  the  single  rotation  parameter  to  characterize  the  space  of 
possible  images. 

Optical  imaging  with  solar  illumination  is  a  bistatic  process.  If  the  object  is  inside  a  sphere  as  in 
the  monostatic  case,  the  receiver  (telescope)  is  at  some  coordinate  pair  and  the  transmitter  (sun)  is  at 
another  one.  Due  to  symmetry,  rotations  of  the  sun  are  not  important.  Hence  rather  than  having  three 
coordinates  specifying  the  orientation  of  the  satellite  we  have  five;  three  specifying  the  satellite  and  two 
specifying  the  sun! 

Low-earth-orbit  objects  can  be  successfully  optically  imaged  just  before  sunrise  and  right  after 
sunset.  A  stable  satellite  at  a  given  elevation  would  experience  solar  illumination  from  the  east  in  one 
case  and  from  the  west  in  the  other.  This  would  cause  different  shadowing  effects.  On  a  given  day  or 
week  there  would  be  essentially  two  illumination  angles  for  a  given  satellite  orientation,  but  throughout 
the  year  these  angles  would  be  changing.  In  addition,  reflections  from  the  earth  affect  the  image.  These 
are  somewhat  dependent  on  weather,  season,  etc. 

In  summary,  when  speckle  imaging  or  speckle  post-processing  adaptive  images  of  astronomical 
objects,  all  of  the  data  can  be  potentially  combined  if  the  (known)  rotation  is  removed  from  the  data. 
When  satellites  are  imaged,  the  data  collected  must  be  processed  in  subsets  characterized  by  up  to  five 
parameters.  Two  of  those  parameters  are  known;  the  satellite  orientation  (motion  solution)  may  also  be 
known  through  independent  means  such  as  radar  imaging.  If  it  is  not,  the  problem  of  determining  the 
three-dimensional  model  of  the  satellite  being  imaged  is  made  more  difficult. 
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We  are  trying  to  characterize  a  three-dimensional  object  using  a  number  of  two-dimensional 
images.  The  number  of  such  images  is  greater  in  the  bistatic  imaging  case.  This  can  have  the  effect  of 
reducing  the  amount  of  data  available  for  integration.  If  this  S/N  issue  is  not  a  problem,  we  will  have 
greater  potential  information  content  as  well  as  a  greater  challenge  to  extract  that  information. 

7.3.3  Field  of  View 

We  have  already  discussed  the  problem  of  nonisoplanicity  in  Section  6.3.4.  The  examples  there 
show  that  objects  subtending  more  than  several  arc  seconds  will  suffer  from  nonisoplanatic  effects 
which  tend  to  reduce  resolution. 

It  would  not  be  unusual  in  optical  space  surveillance  to  be  viewing  a  satellite  at  a  range  of  300  km 
which  is  10  m  in  some  dimension.  In  this  case,  the  corresponding  angular  extent  is  6.9”!  This  is  partic¬ 
ularly  a  problem  for  wavefront  compensation,  since  as  we  have  discussed  speckle  processing  suffers  a 
loss  of  S/N  but  still  has  the  potential  of  recovering  through  integrating  more  data. 

7.4  ANALYSIS  OF  SMALL  AMPLITUDE  AND  PHASE  ERRORS 

Since  we  expect  an  adaptive  feedback  system  to  retain  certain  small  errors,  it  is  natural  to  ask  what 
the  effect  is  of  such  errors.  The  simplest  way  to  characterize  the  error  is  to  determine  its  effect  on  the 
system  MTF.  Indeed,  corrections  for  the  MTF  are  normal  in  speckle  processing. 

To  gain  intuitive  appreciation  for  what  is  going  on,  a  simple-minded  simulation  of  turbulence  was 
produced  in  which  the  phase  standard  deviation  could  be  selected.  While  the  rectangular  phase  error 
spectrum  generated  was  not  as  realistic  as  possible,  we  expect  the  effect  of  small  errors  on  the  MTF  to  be 
qualitatively  realistic.  The  results  were  as  shown  in  Figure  7.7.  The  simulation  is  performed  in  one 
dimension,  for  various  values  of  a^.  The  curve  appears  to  be  similar  to  a  constant  plus  a  sin  jt/jr  at  the 
origin. 

Compensated  imaging  corrects  the  phases  in  a  feedback  loop,  and  as  discussed  above  there  are  sev¬ 
eral  reasons  why  it  cannot  correct  them  perfectly;  that  is  why  the  simulation  modeled  phase  errors.  But 
recall  from  our  discussion  of  turbulence  that  there  is  always  some  scintillation  in  the  magnitudes  as  well. 
Roddier  and  Roddier^  published  an  analytical  analysis  using  certain  approximations  of  the  effect  on  the 
OTF  of  uncorrected  amplitude  errors  due  to  scintillation.  The  curves  which  resulted  look  almost  the 
same  as  Figure  7.7! 

To  understand  why  that  is  so,  we  will  return  to  Goodman."^  His  “Gaussian  random-phase  screen” 
analysis  is  equivalent  to  our  phase  error  scenario.  This  screen  introduces  a  specified  random  phase  shift 
in  the  pupil  plane  of  an  imaging  system  without  disturbing  the  amplitude  of  the  wave. 

The  analysis  is  similar  to  our  previous  discussions  about  the  interferometric  view  of  imaging.  The 
OTF  of  the  screen  is  the  spatial  autocorrelation  function  of  its  amplitude  transmittance  function.  The 
results  shown  by  Goodman  are  again  very  reminiscent  of  Figure  7.7.  We  can  easily  understand  how 
these  similarities  arise  by  recourse  to  communications  theory. 

In  the  theory  of  modulation,  amplitude  modulation  with  carrier  (AM)  and  narrowband  phase  mod¬ 
ulation  (PM)  are  two  well-analyzed  techniques  discussed  in  any  elementary  textbook,  such  as  Lathi.^^  In 
the  communications  context,  a  sinusoidal  carrier  at  frequency  co^  is  modulated  by  an  information  signal 
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Figure  7.7.  Simulated  MTF  of  random  phase  errors. 

f  (/)  for  the  purpose  of  transmitting  /  (r)  to  a  receiver.  In  the  case  of  AM,  the  transmitted  signal  is 

f,it)  =  [A  +  f  (f)]  cos  (Off  (7.2) 

and  in  the  case  of  PM, 

/,(f)  =  /lcos[(Off +  i/(f)].  (7.3) 

The  constant  k  is  known  as  the  modulation  index. 

To  couch  our  optical  problem  in  the  language  of  modulation  theory,  let  the  amplitude  A  be  unity, 
let  (Of  =  0,  and  replace  the  time  function  f(t)  by  some  spatial  error  function  f(x)  of  either  amplitude  or 
phase.  We  then  model  the  spatial  signal  with  amplitude  errors  as 

1+/(A-)  (7.4) 

and  the  signal  with  phase  errors  as 

cos[/(j:)].  (7.5) 

We  are  interested  in  the  autocorrelation  functions  of  these  expressions.  The  simplest  way  to  obtain 
them  here  is  to  first  compute  the  power  spectrum,  and  then  Fourier  transform  (FT).  We  can  appeal  to 
well-known  results  of  modulation  theory  rather  than  doing  any  calculations  ourselves. 

The  AM  analysis  is  fairly  trivial;  clearly  the  power  spectrum  of  Equation  (7.4)  consists  of  an 
impulse  at  the  origin  plus  the  power  spectrum  of  the  modulating  signal,  say  a  Gaussian.  Its  FT  is  a  con¬ 
stant  plus  the  FT  of  a  Gaussian,  which  is  also  Gaussian. 

While  generalized  PM  is  rather  difficult  to  analyze,  involving  numerous  Bessel  functions,  the  spe¬ 
cial  case  of  narrowband  PM  is  much  easier.  Since  we  assume  small  phase  errors,  we  satisfy  the  nar¬ 
rowband  restriction  that  k  f  nil.  It  turns  out  that  the  FT  of  a  narrowband  PM  signal  is  similar  to 
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the  FT  of  an  AM  signal.  The  only  difference  is  that  the  sideband  components  are  multiplied  by  a  scale 
factor  of  jk,  which  amounts  to  a  phase  shift  of  n/2  and  a  scaling  by  the  modulation  index. 

We  conclude  the  power  spectrum  is  the  same  as  in  the  AM  case,  apart  from  scaling,  so  we  have 
explained  the  qualitative  similarity  in  the  OTF  due  to  small  phase  errors  and  small  amplitude  errors.  We 
also  have  obtained  additional  insight  into  the  exact  form  of  the  OTF  of  small  errors.  If  we  can  determine 
the  power  spectrum  of  those  errors,'  the  FT  of  that  spectrum  added  to  a  constant  is  the  desired  OTF. 

Looking  back  at  the  three  examples  which  inspired  this  analysis,  our  simulation  has  a  rectangular 
power  spectrum  which  generated  a  sin  x/x  autocorrelation  function.  This  accounts  for  the  ripples  in  Fig¬ 
ure  7.7.  The  scintillation  paper^  also  shows  ripples  which  reflect  the  spectrum  of  the  amplitude  errors. 
Goodman‘S  uses  a  theoretical  Gaussian,  so  his  curves  show  no  ripples. 

7.5  SELECTING  A  PHASE-RETRIEVAL  ALGORITHM 

Post-processing  of  data  based  on  speckle  techniques  boils  down  to  selecting  a  phase-retrieval  algo¬ 
rithm.  We  have  a  number  of  spatial  frequencies  below  the  seeing  limit  (or  below  the  adaptive  correction 
limit)  whose  phase  is  correct.  Higher  spatial  frequencies  which  we  desire  to  retrieve  have  phases  which 
are  essentially  uncorrelated  with  the  corresponding  spatial  frequency  phases  of  the  object.  We  have 
either  a  good  idea  of  the  magnitude  of  these  frequencies,  differential  phase  information,  or  both.  The 
problem  is  to  determine  the  missing  phases  from  the  available  information. 

As  we  have  discussed  previously  in  Section  6,  there  are  two  basic  approaches  to  the  problem, 
phase  retrieval  from  magnitude  such  as  Fienup’s  method,  and  integration  of  differential  phases  such  as 
Knox-Thompson  or  triple  correlation.  We  have  not  provided  much  guidance  previously  about  which 
methods  to  use. 

The  answer  lies  in  the  well-known  axiom^^-^^  of  image  processing  that  the  phase  is  much  more 
important  than  the  magnitude.  Even  with  no  magnitude  information,  or  the  magnitudes  of  another 
image,  an  image  can  be  recognized.  On  the  other  hand,  with  only  magnitude  information  the  image  is 
not  recognizable. 

We  are  really  estimating  both  the  magnitude  and  the  phase  of  an  image  before  we  display  it.  What 
the  importance  of  phase  tells  us  is  that  we  should  probably  spend  less  time  worrying  about  errors  in  the 
magnitude  since  errors  in  the  phase  are  much  more  important. 

If  we  have  copious  views  of  reference  stars  either  with  very  little  angular  separation  from  the 
object  being  imaged  or  as  part  of  a  very  large  data  set,  we  may  feel  that  our  estimate  of  the  magnitudes 
using  Labeyrie’s  method  is  quite  good.  In  that  case,  following  with  Fienup’s  method  is  probably  a  rea¬ 
sonable  option.  Unfortunately,  that  is  probably  not  the  case  in  most  satellite  imaging. 

If  we  have  substantial  uncertainty  in  the  magnitudes  due  to  a  short  data  set  with  little  suitable  refer¬ 
ence  data,  we  should  be  reluctant  to  base  a  Fienup  reconstruction  on  those  magnitudes.  On  the  other 
hand,  doing  our  best  job  with  one  of  the  differential  phase  algorithms  lets  us  essentially  ignore  the  mag¬ 
nitude  problem.  We  expect  that  if  we  get  good  phases,  the  image  will  contain  detectable  detail  in  those 
spatial  frequencies  even  in  the  face  of  poor  magnitude  estimates.  We  can  then  work  on  improving  the 
magnitudes  without  affecting  the  phases. 
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7,6  OPTICAL  DATA  RESOLUTION  EVALUATION 


The  basic  idea  of  optical  data  resolution  evaluation  is  to  determine  the  quality  of  high-resolution 
optical  data.  There  are  several  reasons  why  this  quality  may  be  uncertain.  Since  the  parameters  of 
speckle  imaging  depend  on  the  seeing  conditions  and  the  nature  of  the  object  being  imaged,  and  we  have 
only  limited  knowledge  of  these  conditions,  even  the  best-planned  measurement  program  has  some 
uncertainty.  If  we  are  post-processing  existing  data  taken  under  less  controlled  circumstances  or  with 
inferior  equipment,  the  uncertainty  will  increase  greatly.  This  situation  might  occur  for  example  if  we 
plan  to  post-process  data  taken  by  an  adaptive  optics  system  or  with  a  camera  of  poor  quality. 

Not  all  objects  contain  significant  energy  at  the  higher  spatial  frequencies.  At  some  frequency 
which  depends  on  the  object  spectrum,  the  object  energy  drops  below  the  noise.  Space  objects  may  be 
rotating  and  introducing  an  unknown  blur  which  affects  the  high  frequencies  first.  These  uncertainties 
remind  us  of  our  discussion  of  reference  stars  in  Section  6. 1.2.2,  We  assume  here  that  there  is  no  handy 
reference  star  so  the  data  must  be  evaluated  using  internal  evidence. 

Assume  we  are  evaluating  a  set  of  n  consecutive  speckle  images  taken  at  a  usefully  short 

exposure  time.  We  first  compute  the  DFTs  of  the  images  ln(u,v)  =  4(f)-  We  would  like  to  know  which 
spatial  frequencies  are  present  and  contributing  to  the  existing  image  quality.  We  then  want  to  go  on  to 
determine  which  differential  spatial  frequency  phases  are  present  and  suitable  for  speckle  phase  recov¬ 
ery.  These  two  measures  are  functions  of  the  magnitude  of  vector  spatial  frequency  |f| ,  so  we  will  aver¬ 
age  both  measures  over  rings  of  constant  I  f  |  in  the  spatial  frequency  plane  to  smooth  the  estimates. 

It  might  be  argued  that  the  first  measure  is  sufficient;  to  determine  the  speckle  improvement  poten¬ 
tial  we  could  simply  perform  the  speckle  processing  and  then  re-evaluate  the  data  using  the  first  quality 
measure.  The  reason  for  proposing  a  measure  of  speckle  potential  is  to  avoid  a  possibly  expensive  and 
complicated  phase  recovery  operation.  The  goal  is  to  obtain  a  relatively  simple  and  sensitive  indication 
of  potentially  recoverable  spatial  frequencies  which  could  guide  sp>eckle  phase  recovery  exp>eriments. 
We  can  also  perform  constant  |f|  averaging  which  is  not  applicable  to  image  phase  recovery;  this  may 
improve  sensitivity  and  enable  us  to  obtain  a  better  estimate  of  spatial  frequency  content  as  a  function  of 

ifi. 

We  adopt  the  following  approach  to  the  problem: 

(a)  Since  Fourier  phases  are  more  important  to  image  quality  than  magnitudes,  we 
base  our  image  evaluation  on  DFT  phase  statistics. 

(b)  Speckle  enhancement  potential  is  based  on  the  statistics  of  differential  DFT  phases, 

(c)  Each  image  quality  measure  is  a  one-dimensional  function  of  spatial  frequency  res¬ 
olution  or  potential  resolution. 

(d)  Since  we  do  not  know  the  true  object  phases  or  “ground  truth,”  we  partition  the 
data  into  two  or  more  disjoint  subsets  and  examine  correlations  within  the  data. 

Closely  correlated  phases  are  assumed  to  correspond  to  true  object  phases. 

We  may  make  use  of  magnitude  information  in  computing  these  phase  measures.  For  example, 
using  only  frequencies  whose  magnitudes  exceed  the  ring  median  or  weighting  the  contributions  would 
probably  improve  performance. 
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7.6.1  Basic  Image  Quality  Measure 
7.6.1. 1  Definition 

Our  basic  image  evaluation  function  5^(|f|)  is  straightforward.  We  start  with  the  assumption  that 
if  the  rms  spatial  frequency  phase  error  is  small  (p)erhaps  less  than  2  n  /  10)^  for  some  frequency 
we  are  getting  reasonable  benefit  from  that  spatial  frequency.  Since  we  don’t  know  the  true 
phases,  we  partition  the  data  into  two  subsets  and  consider  some  function  of  the  phase  difference  at  each 
spatial  frequency  as  an  indication  of  the  mean  squared  error  produced  by  averaging  the  total  data  set.  We 
then  average  over  frequencies  for  which  I  /-/;  I  <  A  (rings)  to  obtain  a  one-dimensional  function  of 
I  f I .  If  we  choose  the  average  squared  phase  differences  as  our  phase  difference  function  we  are 
computing 

^  X  [$,(f)-$2(0]^  (7.6) 

I  f-fi  I  <A 


where  is  the  number  of  discrete  frequencies  within  A  of  and  $,(f)  is  the  average  phase  of  set  /  at 
frequency  f. 

A  more  general  approach  is  to  use  the  sample  phase  variance-,  this  allows  us  to  generalize  by 
dividing  the  data  into  more  than  two  sets.  Dropping  the  bars  for  notational  simplicity,  the  squared  phase 
difference  is 

(4), -4)2)2  =4,2 +4,2 -24),  4)2 


while  the  sample  phase  variance  forK  samples  is 


.2  _  izi _ 

■  K-l 


(7.8) 


where  $  here  is  the  mean  of  the  K  samples.  For  the  case  K  =  2, 


2  ‘!>I  ‘!>2  ,  , 

50  =  —  +  —  -  <t>l  <t>2 


(7.9) 


which  is  just  half  the  squared  phase  difference. 

Using  the  sample  variance  approach  with  K  =  2,  we  obtain  the  image  evaluation  function  j*  ( I  f  I ) 


5|(fi)  = 


1 


N 


ft  I  f-fi  I  <A 


$?(f)  $2(0  . 

— - —  +  -  “  (f)  $2  (0 


(7.10) 


For  this  function,  there  is  nothing  to  be  gained  by  averaging  many  image  DFTs  to  obtain  the  $ 
functions.  As  our  previous  notational  simplification  suggests,  simply  using  two  images  separated  by 
perhaps  30  ms  to  1  s  should  be  quite  suitable.  This  will  not  be  the  case  when  we  discuss  the  speckle 
potential  functions  below. 

We  expect  5^(|f|)  to  look  something  like  Figure  7.8.  There  is  small  phase  difference  until  we 
approach  the  seeing  limit,  at  which  time  there  is  rapid  increase  to  some  limit.  We  can  say  something 
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about  that  limit  in  the  case  in  which  the  two  subsets  have  only  a  single  image  each.  If  the  two  images  do 
not  overlap  in  time,  and  we  are  in  the  noise  regime  or  beyond  the  seeing  limit,  we  expect  the  spatial  fre¬ 
quency  phases  to  be  uncorrelated.  If  the  noise  is  Gaussian,  the  phases  will  be  jointly  uniformly  distri¬ 
buted,  and  the  density  function  will  be  constant  within  the  bounds  of  Figure  7.9(a).  Taking  the  expected 
value  of  Equation  (7.9), 


£[4(|f|)]  =  £ 


2 


-  $1 (0  $2 (0 


(7.11) 


=  £[$1  (OM 


since  the  phases  are  uncorrelated.  This  is  the  variance  of  a  uniformly  distributed  random  variable  over 
the  interval  [-7t,n]  which  is  47t‘/12  or  about  3.29.  We  will  defer  further  discussion  of  this  limit  to  the 
next  section. 


Figure  7.8.  Typical  image  evaluation  functions. 


7.6.1.2  Problems  of  Phase  Discontinuities 

When  we  compare  phases  for  purposes  of  computing  the  sample  phase  variance,  we  take  out  27t 
discontinuities  for  obvious  reasons.  We  can  view  this  process  as  conditionally  adding  or  subtracting  27t 
to  the  second  random  variable.  The  joint  density  function  now  is  constant  within  the  boundaries  of  the 
function  depicted  in  Figure  7.9(b).  This  op)eration  evidently  produces  correlated  random  variables.  The 
cross  term  of  Equation  (7.1 1)  now  is  not  zero,  and  we  must  integrate  to  determine  the  expected  value  of 
this  two  sample  variance.  This  straightforward  calculation  produces  7^12  -nor  about  1.79. 

As  mentioned  above,  we  can  successfully  compute  ^ ^  ( I  f  I )  with  just  two  images  in  which  case  the 
above  analysis  is  appropriate.  When  we  consider  more  than  a  single  image  per  data  set,  which  is 
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(a)  Uncorrelated  phases 


(b)  Correlated  phases 


Figure  7.9.  Joint  probability  density  functions  of  (j)]  (f).  <t>2  (0- 


essential  for  the  differential  phase  function  described  below,  things  get  more  complicated.  We  would 
like  to  take  out  the  phase  discontinuities  in  the  averaging  process  for  each  set.  Etetermining  the  optimum 
set  of  2n  additions  and  subtractions  which  minimize  the  variance  of  a  set  of  data  is  a  non-trivial  optimi¬ 
zation  problem.  We  could  handle  it,  but  doing  it  for  each  spatial  frequency  in  a  large  image  might  be 
burdensome.  The  analysis  of  such  a  process  might  also  be  rather  difficult. 

If  we  proceed  using  a  simple  ad  hoc  generalization  of  our  procedure  for  removing  a  single  discon¬ 
tinuity  from  two  samples,  we  probably  obtain  near-optimum  results.  The  analysis  is  also  straight¬ 
forward,  although  very  tedious.  A  sequence  of  two-dimensional  integrations  is  required  to  analyze  the 
effect  of  this  process  for  white  noise.  Simulating  n  random  images  indicates  that  increases  from 
about  4.1  to  5.4  as  the  n  goes  from  2  to  50.  The  distribution  rapidly  approaches  a  Gaussian  as  one  might 
expect. 

Of  course  our  ad  hoc  discontinuity  removal  process  is  not  completely  defensible  when  dealing 
with  a  Gaussian  distribution  rather  than  a  uniform  one.  A  Gaussian  distribution  with  a*  of  5.4  has  a  of 
2.32,  which  means  that  ±7C  occurs  at  about  1.35  a.  Only  18  percent  of  the  density  exceeds  this  value,  and 
the  density  has  fallen  to  40  percent  of  the  uniform  density  level  at  the  edges.  This  is  not  so  far  from  uni¬ 
form  as  to  introduce  major  deviations  from  the  two  image  case;  indeed  simulations  show  that  the  limit¬ 
ing  value  of  ( I  f  I )  remains  fairly  close  to  1 .8  when  multiple  images  are  processed. 

If  we  omit  the  discontinuity  removal  altogether,  then  random  phases  experience  variance  reduction 
proportional  to  the  number  of  phases  averaged.  Hence  the  expected  distribution  becomes  a  Gaussian 
with  smaller  and  smaller  variance.  This  processing  itself  is  not  suitable  for  a  general  purpose  data  evalu¬ 
ation  function  due  to  its  asymptotic  sensitivity  to  the  exact  numbers  of  images  processed.  We  have  the 
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undesirable  behavior  in  which  the  limiting  value  signifying  the  diffraction  limit  or  the  end  of  useful 
information  varies.  There  are  perhaps  ways  to  deal  with  this  difficulty,  but  phase  discontinuities  would 
degrade  performance  so  it  doesn’t  seem  worth  pursuing. 

7.6.1.3  Interpreting  the  Function 

From  this  analysis  it  follows  that  ^^(Ifl)  should  start  very  small  for  small  spatial  frequencies 
below  the  seeing  limit.  It  should  increase  fairly  rapidly  at  the  seeing  limit  to  an  asymptotic  value  around 
1.8  as  shown  in  Figure  7.8.  If  the  data  are  extremely  noisy,  or  the  object  has  few  significant  spatial  fre¬ 
quencies  beyond  some  point,  the  slope  of  the  increase  may  be  more  gradual  and  the  initial  increase  may 
occur  earlier.  Any  values  of  (If I )  significantly  below  1.8  suggest  some  spatial  frequency  phase  infor¬ 
mation,  but  the  contribution  to  image  quality  is  suspect  unless  5^(|f|)  is  below  0.4  or  0.5.  The  image 
degradation  due  to  phase  errors  depends  on  the  object  and  on  the  magnitude  errors  as  well. 

If  averaging  occurs  in  rings  of  constant  A|  f  | ,  the  number  of  spatial  frequencies  averaged  increases 
with  I  f| .  Hence  we  expect  noisier  estimates  at  lower  I  f  | ;  simulation  clearly  shows  this  effect.  Averag¬ 
ing  in  different-sized  rings  could  compensate  for  this  effect  and  produce  smoother  curves,  although  as 
long  as  the  phenomenon  is  understood  it  is  unlikely  to  cause  much  confusion. 

7.6.2  Speckle  Potential  Evaluation 

7.6.2.1  Phase  Difference  Approach 

The  5^(1  f| )  function  is  relatively  easy  to  compute  and  analyze;  the  phases  are  never  extremely  far 
from  a  uniform  distribution  and  the  lower  frequency  phases  are  very  well  correlated.  We  next  tackle  the 
more  difficult  problem  of  detecting  speckle-recoverable  phases. 

The  idea  behind  the  speckle  potential  evaluation  is  to  compute  phase  differences  with  some  small 
offset  in  the  frequency  plane  in  the  same  way  some  of  the  phase  recovery  algorithms  work.  Rather  than 
doing  anything  with  the  phase  differences,  we  just  statistically  study  them  to  determine  the  potential  for 
phase  recovery. 

The  most  obvious  approach  is  very  similar  to  that  for  the  5^(|f|)  function,  with  DFTs  of  short- 
exposure  images  furnishing  the  raw  data.  However,  instead  of  using  the  phases  of  the  spatial  frequencies 
we  first  compute  the  differences  of  spatial  frequency  phases  separated  by  one  frequency  step  in  u  and  one 
in  V.  The  remaining  processing  follows  the  same  way  as  before.  We  average  a  number  of  frames  to  form 
two  or  more  sets  of  average  differential  phase.  We  then  compute  the  sample  variance  for  each  spatial 
frequency. 

^i$(fi)  =  Tr-  ^  [5$,(f)-A$2(f)]^  (7.12) 

In  Figure  7.10  we  show  the  functions  5^(lf|)  and  calculated  using  two  sets  of  10  simu¬ 

lated  monochromatic  speckle  images  of  a  point  source.  The  MTF  is  simulated  to  roll  off  linearly  to  zero 
at  the  diffraction  limit,  at  which  point  only  noise  contributes  to  the  phases.  In  the  limit  of  completely 
uncorrelated  signals,  5i<j,(|f|)  behaves  the  same  as  5^(|f|),  and  the  asymptotic  value  of  about  1.65  (due 
to  the  multiple  images)  is  approached.  All  of  our  discussion  of  phase  discontinuity  removal,  both  in 
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averaging  the  differential  phases  and  in  computing  the  variance  itself,  applies.  We  see  that  (I f I )  indi¬ 
cates  speckle  phase  recovery  potential  quite  near  to  the  diffraction  limit.  We  see  greater  fluctuations  in 
the  5^  ( I  f  I )  curve  for  smaller  I  f  |  as  expected. 


Figure  7.10.  Evaluation  functions  for  simulated  data. 


This  works  very  nicely  for  simulated  point  sources  but  fails  for  realistic  extended  objects.  The 
problem  occurs  with  the  spatial  frequencies  beyond  the  seeing  limit  and  below  the  diffraction  limit  in 
which  we  have  the  most  interest.  The  analysis  producing  the  asymptotic  value  1 .79  assumes  that  the 
data  are  roughly  uniformly  distributed  over  [-71,71]  and  any  phase  discontinuities  are  removed.  In  the 
case  of  5^  ( I  f  I ),  spatial  frequencies  beyond  the  seeing  limit  are  indeed  independent  and  uniformly  distri¬ 
buted  or  close  to  it.  However,  this  is  not  the  case  for  the  differential  phase  used  to  compute  5^  (I  f| ). 

To  see  this  we  can  consider  a  simple  model  in  which  the  object  is  a  random  field  inside  some 
defined  region  of  the  (x,y)  plane  and  zero  outside.  Goodman  and  Belsher^^  define  a  similar  model  in 
which  the  support  region  is  a  square. 

If  we  have  a  sample  function  of  a  random  field,  and  we  window  it  by  some  support  function,  that  is 
equivalent  to  convolving  its  transform  by  the  transform  of  the  window,  say  a  sine  function  for  a 
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rectangular  window.  The  result  is  that  significant  correlations  appear  between  neighboring  frequencies 
in  the  DFT.  It  follows  that  differences  of  closely  spaced  frequencies  are  smaller  than  they  would  be 
otherwise.  The  differences  are  smaller  in  both  the  real  and  the  imaginary  parts  of  the  complex  FT,  and 
hence  the  phase  differences  are  also  smaller. 

We  now  have  a  phase  distribution  which  is  much  narrower  than  a  uniform  on  [-71, 7i],  so  when  we 
compute  the  variance  for  ( I  f  I )  we  will  automatically  obtain  smaller  numbers.  This  is  independent  of 
whether  or  not  the  spatial  frequency  phase  differences  are  correlated  across  the  two  data  sets.  Clearly  the 
expected  value  of  in  this  region  depends  on  the  support  region  of  the  object  as  well  as  on  the 

speckle  processing  potential.  A  larger  support  region  will  produce  a  smaller  transform  width  and  hence 
produce  less  frequency  domain  correlation. 

We  are  saying  that  the  phase  function  of  a  DFT  is  smooth  due  to  the  support  constraints.  This  does 
not  affect  ^^(Ifl)  since  the  smooth  variation  is  random  from  data  set  to  data  set.  It  does,  however,  pre¬ 
vent  us  from  easily  detecting  speckle-recoverable  frequencies  by  just  looking  for  frequencies  whose  dif¬ 
ferential  phase  variance  is  significantly  below  the  asymptotic  value. 


7.6.2.2  Sample  Correlation  Approach 


An  alternate  approach  is  based  on  sample  correlation  coefficients.  To  avoid  sensitivity  to  the  par¬ 
ticular  distributions,  a  sample  correlation  coefficient  may  be  computed  which  will  always  lie  between  -1 
and  1.  We  can  select  all  of  the  frequencies  in  an  averaging  ring  as  samples;  if  they  show  phase  correla¬ 
tion  between  the  two  data  sets  we  expect  a  positive  correlation  coefficient.  Incidentally,  we  are  again 
restricted  to  partitioning  our  data  into  only  two  sets. 

The  use  of  correlation  coefficients  to  estimate  unknown  signal  and  noise  powers  has  been 
analyzed^^  previously  for  ordinary  one-dimensional  signals.  We  have  a  somewhat  unconventional 
application  for  sample  correlation  coefficients.  Normally  we  would  expect  to  be  given  n  sample  pairs 
from  two,  possibly  correlated,  distributions.  The  number  which  best  represents  the  strength  of  the  linear 
relationship  between  the  two  distributions  is  given  by 

'Z{Xi-x)(yi-y) 

r  =  — -  (7.13) 

ns^Sy 

where  x  and  y  are  the  sample  means,  and  and  Sy  are  the  sample  standard  deviations. 

We  are  proposing  to  apply  Equation  (7.13)  to  a  set  of  sample  pairs  representing  a  number  of  dis¬ 
crete  frequencies  f,-  whose  vector  frequency  magnitudes  are  within  some  A  of  spatial  frequency  magni¬ 
tude  f^  .  The  correlation  we  are  looking  for  is  between  the  data  subsets.  Omitting  the  limits  on  the  sums 
and  the  explicit  4»  dependence  on  f  and  expanding  we  obtain  the  following  function 


(7.14) 


in  a  form  suitable  for  machine  evaluation.  Since  we  know  that  the  phase  means  should  be  zero,  we  can 
explicitly  include  that  knowledge  to  obtain  the  final  expression 
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do  (fi)  = 


(7.15) 


n 


In  cases  where  only  a  few  phase  samples  are  included  in  the  sums,  this  provides  a  substantially  more 
robust  estimator  of  correlation. 

Figure  7.1 1  shows  r^(^(|f|)  calculated  for  the  same  data  used  in  Figure  7.10  except  that  the  point 
source  is  shifted  away  from  the  origin.  (In  the  unique  case  of  a  point  source  at  the  origin,  all  of  the 
object  phases  and  phase  differences  are  zero,  so  the  phase  correlations  are  all  due  to  noise.  Such  shifts 
do  not  affect  ( I  f I ).)  We  can  clearly  extract  the  same  information  from  this  function  that  we  previ¬ 

ously  extracted  from  (I f I )  in  Figure  7.10.  For  more  realistic  extended  objects  (I  f I )  will  retain  its 
general  appearance  while  (I f I )  will  be  smaller,  complicating  the  analysis  as  discussed  above. 


Figure  7.11 .  Sample  correlation  function  for  the  data  of  Fig.  7.10. 


The  asymptotic  correlation  value  is  of  course  due  to  discontinuity  removal.  Using  the  same  type  of 
calculations  as  before,  computing  the  expected  value  of  r\^  (I  f  I )  by  integrating  the  distribution  of  Figure 
7.9(b)  produces  the  value  3/V?7i  or  about  0.68  for  this  limit.  Since  we  use  10  images  per  set,  the  actual 
asymptotic  value  slightly  differs  from  0.68. 

We  are  processing  these  values  as  if  they  were  part  of  a  single  distribution.  The  theoretical 
justification  for  this  approach  is  that  although  we  have  different  distributions  for  different  spatial  fre¬ 
quencies,  we  do  expect  roughly  constant  correlation  between  spatial  frequencies  with  large  magnitudes 
in  a  given  spatial  frequency  magnitude  range.  Those  object  spatial  frequencies  which  happen  to  be  small 
will  have  lower  S/N  and  consequently  lower  sample  correlation.  This  is  why  we  might  want  to  weight 
the  calculation  based  on  DFT  magnitudes  or  delete  the  phases  associated  with  small  magnitudes  from 
the  calculation. 
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7.6.3  Summary 

The  application  of  high-resolution  optics  technologies  to  space  surveillance  is  more  difficult  than 
applying  these  methods  to  astronomy.  In  the  face  of  uncertainty  about  objects  being  imaged,  their 
motion,  and  average  atmospheric  parameters,  data  enhancement  based  on  internal  evidence  becomes 
more  important.  To  guide  this  process,  data  evaluation  methods  based  on  the  same  type  of  internal  mea¬ 
sures  seem  useful. 

The  evaluation  functions  defined  above  represent  a  useful  tool  in  this  process.  Certainly  more  the¬ 
oretical  and  experimental  work  is  needed  to  refine  the  measures  and  their  interpretation.  Additional 
measures  based  on  the  same  general  philosophy  may  also  be  developed  in  the  future.  For  example,  the 
analytical  S/N  calculations  of  Ayers  et  could  be  adapted  to  this  application. 
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8.  SUMMARY  AND  CONCLUSIONS 


We  have  traveled  a  long  road  in  our  effort  to  understand  the  new  technologies  which  support  high- 
resolution  optical  surveillance  of  space  objects.  Starting  from  the  basic  physics  of  solar-illuminated 
imaging  theory  we  have  outlined  the  development  of  these  methods  in  a  quasi-historical  fashion.  The 
main  goal  has  been  to  clearly  explain  the  basic  concepts  with  consistent  notation  while  retaining  compat¬ 
ibility  with  the  references. 

The  reference  list  is  fairly  long  and  reasonably  up-to-date  although  by  no  means  exhaustive.  These 
references  along  with  their  references  constitute  a  more  complete  list  of  relevant  papers.  It  is  certainly 
essential  to  consult  the  literature  to  obtain  a  deef)er  understanding  of  these  topics. 

The  number  of  recent  papers  cited  illustrates  the  great  current  research  interest  in  speckle  imaging. 
It  seems  that  several  very  important  and  interesting  papers  come  out  each  year.  These  are  primarily  from 
the  astronomical  community.  It  is  perhaps  unfortunate  that  the  adaptive  optics  literature  seems  some¬ 
what  disjoint  from  the  speckle  literature;  that  is  gradually  changing  as  adaptive  astronomical  projects  are 
proposed  and  the  commonality  of  the  two  areas  is  realized. 

It  seems  wise  to  pursue  experimental  work  in  space  surveillance  in  both  adaptive  and  speckle  tech¬ 
niques.  The  greatest  gains  will  probably  come  from  the  combination  of  these  two  techniques,  but  useful 
experience  imaging  low-orbit  satellites  can  be  obtained  using  either  method. 
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APPENDIX  A  -  ANALYTIC  SIGNALS 


In  elementary  circuit  analysis  courses,  use  is  frequently  made  of  the  “phasor”  for  linear  systems 
analysis.  A  real  signal  is  represented  as  the  real  part  of  a  complex  signal 

cos  cor  =  Re[e-'“']  (A.l) 


in  order  to  simplify  the  calculations.  If  the  single-sided  spectrum  of  the  complex  phasor  passes  through 
a  linear  system,  taking  the  real  part  of  the  result  produces  the  same  answer  as  working  with  the  two- 
sided  spectrum  of  the  sinusoids. 

This  approach  is  applicable  to  constant  frequencies  only,  so  it  is  natural  to  look  for  a  generalization 
of  the  phasor  which  works  on  arbitrary  waveforms.  Such  a  generalization  exists,  and  is  called  the  ana¬ 
lytic  signal  representation.  This  representation  was  originally  applied  to  the  study  of  modulation,  but 
has  been  found  useful  in  many  other  areas.  Analytic  signals  are  in  fact  analytic  functions  as  defined  in 
complex  analysis,  but  we  need  not  discuss  analyticity  in  that  context.  For  a  more  detailed  look  at  ana¬ 
lytic  signals,  consult  references  such  as  Dugundji,^^  Bedrosian,^^  Van  Trees^’  or  Goodman.^ 

Following  the  idea  of  the  phasor,  the  analytic  representation  \jt(r)  of  the  real  signal  s  (t)  has  the  fol¬ 
lowing  properties: 


s(t)=  Re[\|/(r)] 


W)  = 


/>0 

/<0 


(A.2) 

(A.3) 


where  S(f)  is  the  Fourier  transform  of  sir)  and  'F(/)  is  the  Fourier  transform  of  \|/(r).  The  factor  of  two 
keeps  the  energy  the  same. 

We  can  mention  some  advantages  of  analytic  signals  at  this  point.  First,  due  to  the  one-sided  spec¬ 
trum,  the  convolution  of  one  analytic  signal  with  the  complex  conjugate  of  another  is  zero;  this  greatly 
simplifies  certain  manipulations.  Second,  we  can  immediately  associate  the  envelope  and  phase  of  an 
analytic  signal  with  its  magnitude  and  phase.  These  quantities  are  both  fundamental  in  the  study  of 
modulation;  in  our  optical  application  the  magnitude  is  important. 

At  this  point  we  know  the  real  part  of  the  analytic  representation,  which  is  just  the  signal  itself  We 
need  to  determine  the  imaginary  component.  In  the  frequency  domain, 

W)=  2S(f)uif)  (A.4) 


where  u(f)  is  the  unit  step  function 


/>0 

/<0. 


(A.5) 


Expressing  this  as  a  convolution  in  the  time  domain 


\\i(t)  =  2s(t)* 


S(-0 

2 


1 

jlnt 


(A.6) 


(A.7) 
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The  imaginary  part  of  this  expression  is  known  as  the  Hilbert  transform  of  s{t),  sir).  Hence  the  analytic 
signal  representation  of  s  (t)  is  just  5  (r)  +  jHO- 


We  define  the  complex  envelope  of  the  analytic  signal  \|i(r)  to  be  I  \|i(r)  |  ,  or  s^ ir)+P (r) .  Simi- 

-1  5(0 


larly,  the  instantaneous  phase  of  \|/(/)  is  defined  to  be  tan” 


sir) 


A 
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APPENDIX  B  -  STRUCTURE  FUNCTIONS 


Structure  functions  were  originally  introduced  by  Kolmogorov.  A  couple  of  useful  references  are 
the  works  by  Yaglom^^  and  Tatarskii  both  translated  by  Richard  Silverman.  These  functions  provide  a 
means  of  bringing  the  machinery  develop)ed  for  the  study  of  stationary  random  processes  to  bear  on  cer¬ 
tain  types  of  non-stationary  random  processes.  We  will  not  be  rigorous  in  distinguishing  the  various 
classes  of  stationarity  in  this  Appendix. 

To  motivate  our  discussion,  consider  the  one-dimensional  random  walk.  A  particle  starts  at  the 
origin  and  then  takes  a  random  jump  forward  or  backward  one  unit  each  time  increment.  If  /  (k)  is  the 
position  of  the  particle  at  time  k,  although  E[f  (A:)]  is  zero  for  all  k,  recall  that  the  expected  distance  of 
the  particle  from  the  origin  after  n  steps  is  In  other  words,  the  energy  of  the  process  f(k)  increases 
with  k  so  the  process  is  not  stationary. 

If  we  consider  instead  the  increment  process 

FT{k)^f{k-^T)-f{k)  (B.l) 

we  see  that  this  process  is  stalionar>\  The  random  jumps  exp)ected  in  any  T-step  subintervals  are  statisti¬ 
cally  identical. 

In  this  simple  example  the  process  F  is  exactly  stationary.  In  practical  cases,  such  as  meteorologi¬ 
cal  quantities  in  the  atmosphere,  the  increment  process  will  only  be  approximately  stationary.  The 
approximation  is  good  if  no  appreciable  variation  in  the  process  parameters  occurs  over  the  increment 
time  T,  the  continuous  analog  of  the  discrete  time  T  above. 

Let  /  (r)  be  a  continuous  real  random  process.  Without  assuming  anything  about  the  stationarity 
of  the  increment  process 

F,(t)=f(t^x)-f{t)  (B.2) 

we  want  to  consider  the  autocorrelation  function  of  the  process  F 

RFUut2)^F,{tOF^{t2y  .  (B.3) 

An  algebraic  identity  lets  us  write  this  using  Eq.  (B.2) 

RfU  1  2)  =  y  |<l/  (n  +  T)  -/  (^2)]^  + 

-  <l/ai  +T)-/(r2  +T)]2>-<[/'Oi)-/a2)]0|  (B.4) 

All  four  terms  in  this  expression  are  of  the  form 

=  (B.5) 

so  if  this  structure  function  is  known,  we  can  calculate  the  autocorrelation  function  using  Eq.  (B.4). 
This  is  the  reason  we  work  with  structure  functions. 

If  the  autocorrelation  function  of  a  random  process  is  a  function  only  of  the  difference  of  the  argu¬ 
ments,  then  the  process  is  stationary.  From  Eq.  (B.4)  it  follows  directly  that  if  the  structure  function 
depends  only  on  r,  -  tj,  the  increment  process  is  stationary. 
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It  is  quite  reasonable  to  work  with  increments  of  stationary  processes  as  well.  In  this  case,  we  find 
that  we  can  express  the  structure  function  in  terms  of  the  autocorrelation  function 

D/T)  =  2[;?/0)-;?/t)].  (B.6) 

There  is  much  more  we  could  say  about  structure  functions.  When  investigating  unknown  random 
processes,  it  is  sometimes  convenient  to  assume  stationarity  and  look  at  the  power  spectrum.  If  we  con¬ 
sider  instead  the  increment  process,  we  buy  substantial  protection  against  non-stationarity.  Under  most 
conditions,  if  the  process  turns  out  to  be  stationary  we  can  calculate  the  autocorrelation  function  from  the 
structure  function. 

There  is  one  unnecessary  restriction  we  have  imposed  here;  we  have  assumed  that  the  process  is 
real.  For  remarks  on  the  complex  case,  consult  Yaglom.^^ 
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